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APPLICATIONS  OF  THE  HOMOMORPHIC  TRANS F0R21AT ION 
TO  TE'IE  DOMAIN  MEASUREMENT  PROBLEMS 


Sedkl  M.  Rlad  and  Norris  S.  Nahman 


This  report  presents  a  study  of  the  theory  and  application  of  the  ho~o- 
morphic  transformation  to  deconvolution  problems  occurring  in  time  domain 
measurements  in  the  picosecond  time  domain.     A  homomorphic  deconvolution  trans- 
form was  developed  and  applied  successfully  to  remove  the  time-windowing  restric- 
tion required  in  many  time  domain  measurements.     Examples  were  presented 
including  problems  in  time  domain  analysis  of  linear  networks  and  dielectric 
spectroscopy,  and  scattering  and  multiple  reflection  in  antenna  (radiation) 
systems  were  considered  and  treated.     Also  considered  was  the  development  of  a 
model  for  a  28  picosecond  resolution  feed-through  sampling-head,  and  the 
model's  step  response  was  computed.     Simulation  studies  were  performed  using 
typical  input  waveforms  and  the  oscilloscope  model.     The  homomorphic  transfor- 
mation was  used  to  deconvolve  the  model's  impulse  response  from  the  simulated 
output.     Comparison  of  the  deconvolved  output  waveforms  with  the  input  waveforms 
showed  agreement  within  the  accuracy  of  the  sampled-data  simulation. 

Key  words:  Antenna  scattering;  deconvolution;  dielectric  spectroscopy; 
homomorphic  transformation;  modeling  of  sampling  gates;  signal  process- 
ing; time  domain  measurements;   time  domain  ref lectometry . 

1 .  INTRODUCTION 

This  report  is  concerned  with  the  study  of  the  theory  and  application  of  the 
homomorphic  transformation  to  deconvolution  problems  occurring  in  time  domain  spectroscopy 
and  scattering.     Typical  time  domain  problems  that  fall  under  this  category  include  the 
following:     dielectric  and  magnetic  spectroscopy,  time  domain  ref lectometry  and  transmis- 
sion (network  spectroscopy),  and  wave  scattering  and  multiple  reflections  in  antenna 
problems.     Conventionally,  time  domain  measurements  methods  are  based  upon  and  limited  by 
time-windowing  to  select  the  relevant  or  desired  waveform  epoch. 

In  this  report  a  new  method  is  developed  based  upon  the  mathematical  homomorphic 
transformation  which  removes  the  time-windowing  restriction  required  in  the  usual  time 
domain  measurement  method.     This  allows  the  use  of  overlapping  time  vrindows  which  effec- 
tively removes  the  time-window  restriction  and  thus  engenders  a  new  method  for  time 
domain  spectroscopy,  and  in  general,  a  new  method  for  time  domain  measurements  in  the 
picosecond  time  domain  [1] . 

Furthermore,  because  time  domain  representations  have  the  property  of  directly 
exhibiting  reflected  or  scattered  contributions,   this  method  has  significant  impact  in  the 
broad  class  of  problems  falling  under  the  category  of  scattering  problems  which  are 
encountered  widely  in  electrical  engineering  and  in  physics   [1] . 

Specifically,  this  report  develops  four  problems,  or  applications,  which  are  encountered 
in  various  forms  in  time  domain  spectroscopy: 


1.       Deconvolution  of  incident  and  scattered  components  in  the  time  domain 
analysis  of  linear  networks . 
-  1     ,   2.      Removal  of  the  time  windowing  restriction  caused  by  a  finite-sized  specimen 
in  time  domain  spectroscopic  studies  of  dielectrics. 

3.  Removal  of  scattering  components  from  antenna  test  range  data. 

4.  Determination  (deconvolution)  of  the  excitation  waveform  from  the  observed 
"  sampled-data  response  of  a  picosecond  domain  sampling  oscilloscope. 

In  summary,  this  report  explores  in  depth  the  theory  and  application  of  the  homo- 
morphic  transformation  to  scattering  problems  using  the  above  topics  as  examples,  including 
experimental  and  computer  data  reduction  studies. 

This  report  is  organized  into  seven  chapters  and  three  appendices.     The  present 
introduction  being  the  first  chapter. 

Chapter  2  will  be  devoted  to  the  study  of  the  theoretical  aspects  of  the  homomorphic 
transformations,  and  in  particular,  the  homomorphic  deconvolution  transformation.    Also  the 
computational  realization  of  the  homomorphic  deconvolution  transformation  will  be  given. 

In  chapter  3  the  homomorphic  deconvolution  transform  developed  in  chapter  2  will  be 
applied  to  various  time  domain  linear  network  analysis  problems.     The  applications  will 
cover  both  time  domain  ref lectometry  and  transmission  measurements  for  pure  resistive 
networks  (e.g.,  commensurate  transmission  lines  network),  as  well  as  networks  with  general 
(reactive)  reflective  elements  (discontinuities) . 

Chapter  4  will  present  the  theory  and  limitations  concerning  time  domain  dielectric 
spectroscopy.     The  homomorphic  deconvolution  transformation  will  be  applied  to  remove  the 
restrictions  of  the  conventional  method. 

In  chapter  5  the  homomorphic  deconvolution  transform  will  be  applied  to  another  class 
of  time  domain  scattering  problems:     the  removal  of  scattering  components  from  antenna 
test  range  data. 

Chapter  6  will  be  concerned  with  the  transfer  network  identification  (modeling)  of 
feed-through  sampling  head.     The  theory  and  construction  of  the  sampling-head  will  be 
reviewed.     Various  experimental  and  computational  steps  involved  in  the  modeling  process 
will  be  presented.     The  model  will  be  used  to  compute  the  sampling-head  step  response.  A 
deconvolution  routine  is  to  be  developed  in  order  to  deconvolve  the  sampling-head's 
response  contribution  from  the  observed  waveform  on  the  oscilloscope  to  get  an  estimate  of 
the  actual  waveform  at  the  input  port  of  the  oscilloscope. 

Chapter  7  will  contain  the  summary  and  conclusions  of  the  work  presented  in  this  report. 

Relevant  detailed  mathematical  developments  and  computer  programs  are  referred  to  in 
the  text  and  presented  in  appendices  A  through  C. 
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2.     THE  HOMOMORPHIC  TRANSFORMATION 


2.1  Introduction 

It  is  often  encountered  in  time  domain  measurements,  as  well  as  in  other  various 
applications,  that  an  observed  signal  or  a  waveform  x(t)  can  be  represented  by  the  convolu- 
tion of  two  (or  more)  components  x^(t)  and  X2(t);  that  is, 

GO 

x(t)  =  /    x^(t')  X2(t-t')  dt'  (2.1a) 

—CD 

or  using  the  simplified  notation 

x(t)  =  x^(t)  *  X2(t)  (2.1b) 

where  the  (*)  denotes  the  convolution  process.     In  the  frequency  domain,  eq.    (2.1)  takes 
the  form 

X(ej'')  =  X^(ej")  •  X2(e>)  (2.2) 

where  X(e-^"),  X^(e-''^),  and  X2(e-''^)  are  the  Fourier  transforms  of  x(t),  x^(t),  and  X2(t), 
respectively,  all  complex  functions  of  e-''^. 

Some  of  the  examples  of  time  domain  convolutions  are  (a)  the  reflected  signal  on  a 
transmission  line  which  is  the  convolution  of  the  incident  signal  and  the  impulse  response 
of  the  reflecting  discontinuity  and  (b)  the  observed  signal  using  an  oscilloscope  which  is 
the  convolution  of  the  measured  signal  (at  the  input  gate)  and  the  oscilloscope's  impulse 
response . 

In  order  to  separate  the  components  of  a  convolution,  a  deconvolution  (separation) 
process  is  needed.     In  the  past  several  deconvolution  methods  were  developed  for  various 
convolution  problems;  they  all  required  the  knowledge  of  either  one  of  the  convolution 
components,  and  in  some  cases,  information  about  the  signals  being  analyzed.     A  fairly  new 
method  of  signal  processing,  the  homomorphic  transformation  [2,3,4,5],  will  be  described 
in  the  present  chapter  and  will  be  used  throughout  the  following  chapters  in  deconvolving 
time  domain  signals.     The  new  method  does  not  require  the  knowledge  of  any  of  the  convolu- 
tion components  and  suits  the  class  of  problems  where  the  frequency  domain  forms  of  the 
convolution  components  have  substantially  different  rates  of  variation.     As  will  be  shown 
later,  the  time  domain  ref lectometry  and  spectroscopy  problems  do  yield  to  the  new 
deconvolution  method. 

Since  all  signal-processings  that  will  be  encountered  in  the  deconvolution  process 
will  be  carried  out  using  digital  computers,  the  analysis  will  be  presented  in  the  discrete 
(sampled-data)  form. 
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■    2.2    The  Homomorphic  Transformation 

A  homomorphic  system  H,  figure  2.1,  is  a  nonlinear  system,  which  obeys  the  generalized 
principle  of  superposition  [2,3,4];  i.e., 

H[x^(n)  0^  x^Cn)]  =  H[x^(n)]  0^  H[x2(n)]  (2.3) 

and 

H[c  :  x^(n)]  =  c  i  H[x^(n)]  (2.4) 

where  x^(n)  and  X2(n)  are  two  input  sequences  to  the  system  H,  n  is  the  sequence  variable, 
c  is  any  scalar  quantity,  0^  (or  0^)  denotes  a  rule  for  combining  inputs  (or  outputs)  with 
each  other  (e.g.,  addition,  multiplication,  convolution,  etc.),  and  :  (or  \)  denotes  a 
rule  for  combining  input  (or  output)  with  a  scalar.     The  class  of  linear  systems  is  a 
special  case  for  which  0^  and  0^  are  addition  and  :  and  i  are  multiplication.  Such 
systems  are  said  to  have  an  input  operation  0^  and  an  output  one  0^. 

The  canonic  representation  of  homomorphic  systems  is  shown  in  figure  2.2.  Both 
2r 

(2.4), 


systems        ,  and  D^^  obey  the  generalized  principle  of  superposition,  eqs .   (2.3)  and 
i  o 


The  input  operation  for         is  0^  while  the  output  one  is  +  (addition) ;  the  system  L 

i  -1 

is  a  linear  system,  and  the  system  D_     transforms  from  +  to  0  . 

•^  0  o 

o 

'  -  ^         ■  -  oi  ;      :      IV  -  ;  :  ;  2.3    The  Separation  Idea 

Referring  to  figure  2.2  and  letting  x(n)  denote  the  combination  of  the  input  sequence 
Xj^(n)  and  X2(n),  it  can  be  written  that 

[x(n)]  =  Dq  [x^(n)  0.  X2(n)]  =        [x^(n)]  +  [x2(n)] 
i  i  i  i 

'     '  =  x^(n)  +  X2(n) 


=  x(n)  (2.5) 


L[x(n)]  =  L[xj^(n)  +  X2(n)]  =  L[x^(n)]  +  L[x2(n)] 


=  y(n)  (2.6) 


DQ^[y(n)]  =  DQ^[yi(n)  +  y2(n)]  =  DQ^[y3^(n)]  0^  DQ^[y2(n)] 
o  o  o  o 


=  y^Cn)  0^  y2(n) 

=  y(n)  .  (2.7) 


Now,  if  it  is  required  to  recover  one  of  the  input  signals  (e.g.,  x^(n))  from  the  combined 
sequence  x(n),  the  linear  system  L  must  have  the  filtering  property 

L[x^(n)  +  x^Cn)]  =  x-j^(n);  (2.8a) 

i.e. , 

y(n)  =  x^(n).  (2.8b) 
Then,  choosing 

Dq    =  %  (2.9) 
o  1 

eqs.   (2.7)  and  (2.5)  yield 

y(n)  =  DQ-'-[Xj^(n)]  =  x^(n) .  (2.10) 
i 

In  conclusion,  the  separation  of  the  combined  sequences  x^(n)  and  X2(n)  depends  on  the 
success  in  separating  their  transforms  x^(n)  and  X2(n) ,  which  depends  in  turn  on  the 
properties  of  the  signals  involved. 


2.4    The  Homomorphic  Deconvolution 

As  mentioned  earlier  in  section  2.1,  the  deconvolution  is  the  process  of  separating 
convoluted  signals  from  each  other.     The  homomorphic  deconvolution  is  then  the  use  of  a 
suitable  homomorphic  system  (of  transformations)  to  separate  (deconvolve)  a  convolution. 
The  input  and  output  operations  of  a  homomorphic  deconvolution  system  are  the  convolution 
operation  (*) ,  according  to  eq.   (2.9).     The  canonic  form  of  a  homomorphic  deconvolution 
system  is  shown  in  figure  2.3.     In  the  figure,  the  transform        converts  the  convolution 
operation  into  an  addition.     The  transform  D^''"  is  the  inverse  transform  of  D^^.     In  con- 
structing the  system        to  do  its  function,  there  are  two  useful  mathematical  operations 
that  will  be  used : 

a.  The  z-transform  z[     ],  appendix  A,  converts  a  convolution  into  a  product, 

eq.   (A. 13).     In  other  words,  the  z-transform  is  a  homomorphic  transform  with  (*) 
and  (•)  as  the  input  and  output  operations,  respectively. 

b.  The  complex  logarithm  operator  log  converts  a  product  into  a  sum;  i.e.,  the  log 
is  a  homomorphic  transform  with  (•)  and  (+)  as  the  input  and  output  operations, 
respectively. 

The  simplest  form  for  the        transform  can  then  be  the  cascade  combination  of  the  z- 
transform  and  the  complex  logarithm,  figure  2.4a,  but  if  it  is  preferable  in  the  separation 
process  to  deal  with  sequences  rather  than  their  z-transf orms ,  the  addition  of  an  inverse 
z-transform  z  ^[     ],  eq.   (A. 3),  which  is  a  linear  one,  will  not  alter  the  output  operation, 
figure  2.4b. 
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Studying  the  signal  flow  in  the  constructed  system  D^^  of  figure  2.4,  with  x(n)  as  the 

convolution  of  two  sequences  Xj^(n)  and  X2(n),  the  following  relations  are  obtained: 

X(z)  =  z[x(n)]  (2.11a) 

=  z[x^(n)  *  x^Cn)]  (2.11b) 

=  z[x^(n)]  •  z[x2(n)]  (2.11c) 

■                  =  X^(z)   .  X^(z)  (2. lid) 

X(z)  =  log[X(z)]  (2.12a) 

=  log[X^(z)   •  X2(x)]  (2.12b) 

=  log[X^(z)]  +  log[X2(z)]           ,  (2.12c) 


X^(z)  +  X2(z) 


(2.12d) 


x(n)  =  z'^[X(z)]  (2.13a) 

=  z"^[X^(z)  +  X^(z)]  (2.13b) 

=  z"^X^(z)]  +  z"^X2(z)]  (2.13c) 

=  x^(n)  +  X2(n)   .  (2.13d) 

The  reason  for  denoting  the  sequence  variable  in  eq.   (2.13)  by  n  is  to  distinguish  the 
new  domain  from  the  original  n-domain.     This  n-domain  is  a  nonphysical  domain;  it  is  a  pure 
mathematical  one  or  what  is  to  be  called  the  hypothetical  n-domain.     If  the  n-domain  is  the 
time  domain,  then  the  n-domain  is  to  be  referred  to  as  the  hypothetical  time  domain.  The 
signal  x(n)  in  the  hypothetical  n-domain  is  historically  named  the  cepstrum  (reversing  the 
order  of  the  first  four  letters  of  the  word  spectrum). 

The  separation  process  will  then  be  achieved  using  linear  filters  in  the  n-domain  to 
split  x(n)  into  its  components  x^(n)  and  X2(n).     The  filtration  problem  will  be  greatly 
simplified  if  x(n)  is  real.     Assuming  that  both  x(n)  and  x(n)  are  real  stable  sequences  is 
not  in  fact  a  restriction  since  these  are  the  properties  of  the  majority  of  signals  encountered 
in  this  work.*    Besides,  simple  mathematical  operations  can  change  the  stability  conditions 
of  a  given  signal.     For  example,  exponential  weighting  can  be  used  to  secure  stability  as 
being  discussed  in  appendix  A.     This  assumption  implies  that  the  regions  of  convergence  of 
both  X(z)  and  X(z)  must  include  the  unit  circle,  see  appendix  A. 

On  the  unit  circle  z  =  e~''^,  X(z)  can  be  expressed  as  X(e-''^) , 

X(ej")  =  X^(e^-)  +  jX^(ej-)  ^^.14) 


*Here  a  stable  signal  is  defined  as  one  whose  z-transform  contains  no  poles  exterior  to 
or  on  the  unit  circle.     Note  that  this  definition  corresponds  to  a  stable  signal  as  being 
defined  as  one  whose  Laplace  or  s-transform  contains  no  poles  on  the  imaginary  axis  or 
in  the  right-hand  plane. 
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where  the  subscripts  R  and  I  refer  to  real  and  imaginary  components  of  a  complex  quantity. 
Since  x(n)  is  assumed  to  be  real,  X_,(e-''^)  and  X  (e-^*^)  must  respectively  be  even  and  odd 
periodic  functions  of  w  with  period  2tt.     Both  functions  must  be  analytic  (continuous 
functions  of  oo)  to  meet  the  stability  requirements.     Since  X  =  log  X  by  virtue  of  eq.  (2.12), 
the  analyticity  of  the  complex  logarithmic  function  is  to  be  studied. 
The  complex  logarithmic  function  log  is  defined  as 

log  X  =  log|x]  +  j  arg[X]  (2.15a) 
-  ir  <  arg[X]  <_  +  n.  (2.15b) 

Using  eqs.  (2.12)  and  (2.14),  together  with  the  above  definition,  the  real  and  imaginary 
components  of  X(e'^'^)  are 

X^(ej'^)  =  log|x(e^'^)|  (2.16a) 

and 

X^(ej")  =  arg[X(e^'')]  .  (2.16b) 

The  real  part,  eq.   (2.16a),  is  analytic  as  long  as  the  z-transform  X(z)  does  not  have  any 
zeros  on  the  unit  circle.     However,  the  exponential  weighting  method,  appendix  A,  can  be 
used  to  move  any  zeros  lying  on  the  unit  circle  to  make  eq.   (2.16a)  analytic.     The  imaginary 
part,  eq .   (2.16b),  is  not  a  continuous  function  according  to  eq.   (2.15b);  this  does  have 
two  undesired  effects: 

a)  The  eq.   (2.12c)  will  not  (always)  be  true.     Using  eq .   (2.14),  eq .   (2.12)  becomes 

log|x|  =  log|X^|  +  loglx^l  (2.17a) 

arg[X]  =  arg[X^]  +  arg[X2]  (2.17b) 

where  the  three  arguments  in  eq.  (2.17b)  are  all  defined  in  the  region  -tt  to  +tt. 
The  relation  of  eq.  (2.17a)  is  always  true,  but  eq .  (2.17b)  is  not  because  of  the 
restrictions  on  the  arguments. 

b)  X(z)  will  not  (always)  be  a  valid  z-transform  because  of  the  noncontinuity  of  the 
imaginary  component  of  X(=  log  X),  thus  invalidating  eq.  (2.13). 

The  solution  to  the  complex  logarithm  problem  is  simple;  it  is  to  replace  the  dis- 
continuous complex  logarithm  by  the  continuous  complex  logarithm  clog  defined  as 

clog[X]  =  log|x|  +  j  carg[X]  (2.18a) 

carg[X]  =  arg[X]  +  2mTr  (2.18b) 

where  carg  is  the  continuous  argument  function.     This  new  function  is  formed  by  recon- 
structing the  argument  function  to  be  a  continuous  one  by  addding  (or  subtracting)  integer 
multiples  of  2-n ,  figure  2.5. 

-7- 


In  conclusion,  for  the  system  D^^  to  function  properly  (i.e.,  convert  the  real  stable 
convolution  x(n)  into  a  real  stable  sum  x(n)),  the  constructed  clog  function  must  be  used 
instead  of  the  log  function;  and  the  real  and  imaginary  components  of  both  X(e-''^)  and 
X(e'^'*^)  must  be  respectively  even  and  odd  continuous  periodic  functions  of  cd  with  period  Zir. 

The  inverse  system  D^"''  should  perform  the  reverse  processes  of  the  system  D^;  i.e., 
convert  an  addition  at  the  input  to  a  convolution  at  the  output.     D^''"  is  constructed  by 
reversing  the  order  of  the  transforms  of        replacing  each  transform  by  its  inverse, 
figure  2.6.     The  system  D^'''  does  not  have  any  ununiqueness  problems.     Since  the  system 
D^'''  is  the  inverse  transform  of  D^, 

D;^D^[x(n)]]  =  x(n).  .      ^  . 

Since  both  x(n)  and  x(n)  are  assumed  to  be  real  and  stable,  it  follows  that  both  y(n) 
and  y(n)  are  both  real  and  stable.     Thus  the  regions  of  convergence  of  both  Y(z)  and  Y(z) 
must  include  the  unit  circle. 


2.5    Computational  Realization  of  the  Homomorphic  Deconvolution  Transform 

In  most  physical  situations  the  observed  signals  have  to  be  dealt  with  numerically 
rather  than  analytically.     The  digital  computer  is  often  used  to  perform  various  kinds  of 
numerical  processing  on  such  signals.     For  the  separation  of  physical  convolution  signals, 
a  digital  computer  version  of  the  homomorphic  deconvoluting  transform  has  to  be  derived. 

In  section  2.4,  it  was  decided  that  the  z-transforms  X(z),  X(z) ,  Y(z),  and  Y(z)  have 
to  include  the  unit  circle  in  their  regions  of  convergence.     Thus  the  inverse  z-transforms 
can  then  be  calculated  using  the  unit  circle  as  the  integration  contour,  eq.   (A. 3).  Ac- 
cordingly, the  evaluation  of  the  above  mentioned  z-transforms  is  only  needed  on  the  unit 
circle;  i.e.,  it  is  employed  to  calculate  the  signal  Fourier  transforms  rather  than  the 
z-transforms,  eq .   (A. 6).     Because  of  the  discreteness  and  finiteness  required  in  numerical 
computational  methods,  the  input  signal  has  to  be  used  in  its  discrete  version.  The 
discrete  Fourier  transform,  figure  2.7,  shows  the  computational  realization  of  the  homomorphic 
deconvoluting  transform,  in  which  the  signals  are  related  as  follows: 


X(k)  =  XCeJ'") 
N-1 


w  =  2Trk/N 
-j  (2TTk/N)n 


=    I    x(n)e-J^"^'^'^^"  (2.19) 
n=0 

X(k)  =  clog[X(k)] 

=  log|x(k)  I  +  j  carg[X(k)],  (2.20a) 
carg[X(k)]  =  arg[X(k)]  +  m2TT.  (2.20b) 
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Here  m  may  be  zero  or  a  positive  or  negative  integer,  such  that  carg[X(k.)]  is  continuous. 

kh  =1  Y  X(k)ej(2"^/^>"  (2.21) 
k=0 

y(n)  =  L[x(n)]  (2.22) 

Y(k)  =  Y^(n)  e-J^'^^/")"  (2.23) 
n=0 

Y(k)  =  exp[Y(k)]  (2.24) 

y(n)  =iT^(k)e^^'"^^''^"  •  (2.25) 
k=0 

It  has  to  be  noted,  that  the  input  sequence  x(n)  is  the  convolution  of  two  sequences  x^(n) 
and  X2(n)  both  of  finite  number  of  samples  N;  consequently,  the  cepstrum  x(n)  is  the  sum  of 
the  two  contributions  x^(n)  and  X2(n)  due  to  x^(n)  and  X2(n),  respectively.     The  output  of 
the  linear  filter  y(n)  should  be  equal  to  either  of  the  components  Xj^(n)  or  X2(n)  to  yield 
a  system  output  y(n)  which  correspondingly  equals  either  Xj^(n)  or  X2(n).     Also,  it  is  to  be 
noted  that  all  the  sequences  involved  in  the  transform  equations,  eqs.   (2.19)  through 
(2.25),  either  in  the  n,  k,  or  n-domains  are  each  of  a  finite  number  of  samples  N,  and  they 
each  satisfy  the  relation 

h(m)  =  h(m  +  rN) ;  r  integer;  (2.26) 

i.e.,  all  the  computations  need  to  be  done  only  for  sequence  variables  ranging  between  0 
and  N-1. 

The  system  of  figure  2.7  was  simulated  by  the  computer  subroutine  HDT  given  in 
appendix  B.     In  cases  where  the  input  and  output  signals  to  the  deconvolution  routine  are 
represented  by  their  discrete  Fourier  transforms,  the  subroutine  HDF  is  to  be  used.  HDF  is 
the  same  as  HDT  except  that  it  does  not  contain  the  first  DFT  and  the  last  DIET. 

2.6    Computational  Hint;  the  Linear  Phase  Elimination 

Let  the  sequence  x(n),  the  input  to  the  homomorphic  deconvolution  transform  of 
figure  2.7,  be  the  result  of  delaying  another  sequence  s(n)  by  n^  samples;  i.e., 

x(n)  =  s(n-n^).  (2.27) 
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Following  the  transform  equations,  eqs .   (2.19),   (2.20),  and  (2.21),  it  can  be  written 


that 


X(k)  =  S(k) 


-j(2Trk/N)n^ 


(2.28) 


e 


X(k)  =  S(k)  -  j(27rk/N)n^  , 


(2.29) 


and 


x(n)  =  s(n)  +  (-1) 


n+1 


n  /n  . 

o 


(2.30) 


The  second  term  in  eq .   (2.30)  is  of  the  form  of  a  decaying  oscillation.     This  term  is 
called  the  linear  phase  term  because  it  is  due  to  the  linear  phase  component  in  the  Fourier 
transforms,  eqs.,   (2.28)  and  (2.29).     Assuming  that  it  is  desired  to  perform  linear  filtering 
on  s(n),  the  linear  phase  term  must  be  removed.     Said  another  way,  s(n)  contains  the  desired 
information  while  the  addition  of  the  linear  phase  term  obscures  this  information.  Ac- 
cordingly, the  linear  phase  term  of  x(n)  must  be  removed.     Figure  2.8  illustrates  how  the 
linear  phase  term  distorts  the  information  s(n)  contained  in  x(n) . 

The  existence  of  such  a  linear  phase  term  also  contradicts  the  requirements  of  having 
odd  periodic  continuous  phase  (argument)  function  with  co.     This  linear  phase  component  has 
to  be  eliminated,  and  the  best  place  to  perform  the  elimination  is  right  after  taking  the 
clog  of  the  DFT  of  the  signal.     Figure  2.9  shows  the  linear  phase  elimination  process. 


This   chapter  started  with  a  review  of  the  convolution  relations  in  both  the  time  and 
frequency  domains,  and  the  need  for  deconvolution  (separation)  techniques. 

The  general  class  of  the  homomorphic  systems  was  defined,  and  their  properties  were 
stated  showing  how  they  can  be  used  for  signal  separations.     In  particular,  the  homomorphic 
deconvolution  transform  was  developed.     The  related  problems,  such  as  the  nonanalyticity  of 
the  complex  logarithmic  function  and  the  effect  of  linear  phase  component,  were  presented 
and  discussed,  and  solutions  to  such  problems  were  given. 

The  computational  realization  of  the  deconvolution  transformation  was  made  possible  in 
the  form  of  digital  computer  subroutines  HDT  and  HDF. 

In  the  next  chapters,  the  new  deconvolution  method  will  be  put  into  use  in  deconvolving 
various  kinds  of  time  domain  signals,  such  as  those  observed  in  time  domain  spectroscopy 
measurements . 


V 


2 . 7  Summary 
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x(n)  =  x^(n)0.X2(jn 


1  

(n) 

H 

y(n)  =  H[x(n)] 
'=  H[x^(n)]0^H[x2(n)] 


Figure  2.1     A  Homomorphic  System  H  with  0^,   and  0^ 
as  the  Input  and  Output  Operations, 
Respectively. 


x(n) 

+  + 

+  + 

0, 

0 

0 

! 

L 

y(n) 

y(n) 


Figure  2.2.    Canonic  Representation  of  Homomorphic 
Systems . 


x(n) 


+  + 

x(n) 

+  + 


y(n) 


y(n; 


Figure  2.3.    The  Canonic  Form  of  a  Homomorphic 
Deconvolution  System. 
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* 

»  • 

+ 

z[  ] 

log 

x(n) 

X(z) 

X(Z) 


a.     The  Output  Signal  is  a  z-Transf orm. 


* 

z[  ] 

log 

+       +  ■ 

x(n) 

X(z) 

X(z) 

b.     The  Output  Signal  is  a  Sequence. 


Figure  2.4.  Possible  Forms  for  Construction  of  the 
Transform       ,  a  Homomorphic  System  with 
(*)   as  the  Input  Operation  and   (+)  as 
-    the  Output  Operation. 
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arg[X(ej*")] 


TT 


M 

 ^ 

V 

7T  W 

 ► 

71  U) 

carg[X(e^'")] 


Figure  2.5.  The  Reconstruction  of  the  Argument 

Function  arg  to  Create  the  Continuous 
carg  Function. 


-13- 


y(n) 


z[  ] 


EXP[  ] 

*  * 

z-\  ]  1 

Y(Z) 

y(n) 


Figure  2.6.  The  Transform  D 
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LINEAR 
FILTER 
L[  ] 


y(n) 


DIFT[  ] 


EXP[  ] 


DFT[  ] 


Figure  2.7.   The  Computation  Realization  of  the 

Homomorphic  Deconvolution  Transform. 

DFT  =  Discrete  Fourier  Transform. 

DIFT  =  Discrete  Inverse  Fourier 
Transform . 

CLOG  ^  Continuous  Complex  Logarithm 
Function. 
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Figure  2.9.    The  Elimination  of  the  Linear  Phase 
Component . 
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3.     THE  THEORY  AND  APPLICATION  OF  THE  HOMOMORPHIC  TRANSFORMATION 
TO  DECONVOLUTION  IN  TIME  DOMAIN  LINEAR  NETWORK  ANALYSIS 


3.1  Introduction 

In  chapter  2,  the  homomorphic  deconvolution  transform  was  developed.    As  it  was  shown, 
the  separation  (deconvolution)  idea  was  based  on  the  property  of  this  transform  to  convert 
a  convolution  of  two  sequences  (waveforms)  into  a  sum  of  their  contributions  in  a  hypo- 
thetical time  domain  called  the  cepstrum.     Separating  filters  are  then  to  be  applied  to  the 
cepstrum  sequence  to  select  only  one  contribution  at  a  time.     The  separation  can  be  realized 
easily  if  the  two  cepstrum  contributions  were  discernible  from  each  other;  i.e.,  each  one 
can  be  viewed  in  a  separate  window,  in  which  case  simple  (ideal)  bandpass  type  (mathematical) 
filters  can  be  used  for  the  separation  process. 

In  the  present  chapter,  a  class  of  problems  that  meet  the  above  mentioned  requirement 
for  easy  separation  will  be  treated;  this  class  of  problems  is  the  analysis  of  linear 
microwave  networks  using  time  domain  techniques. 

■  3.2    Time  Domain  Network  Analysis 

The  time  domain  method  of  network  analysis  is  based  upon  observing  the  reflections  (or 
transmission)  of  an  incident  time  domain  waveform  back  from  (or  through)   the  network  to  be 
analyzed.     The  observed  signal  contains  information  on  both  the  incident  signal  and  the 
networks  impulse  response.     The  typical  experimental  setups  for  both  the  reflection  and 
transmission  cases  are  given  in  figure  3.1.     The  method  of  observing  the  reflected  waveform 
is  commonly  called  the  TDR  (time  domain  ref lectometry)  method;  on  the  other  hand,  the 
abbreviation  TDT  will  be  assigned  to  the  time  domain  transmission  measurement  method.  In 
TDR,  figure  3.1a,  the  observed  reflected  waveform  can  be  viewed  as  the  convolution  of  the 
incident  waveform  and  the  reflection  impulse  response  of  the  network.     While  in  the  TDT, 
figure  3.1b,  the  transmission  waveform  is  the  convolution  of  the  incident  waveform  and  the 
transmission  impulse  response  of  the  network.     The  general  nature  of  the  impulse  responses 
of  the  microwave  networks  that  will  be  considered  in  this  chapter  is  of  the  form  of  a  pulse 
train,  because  of  discretely  distributed  discontinuities  along  the  network  transmission 
line.     The  frequency  domain  form  of  such  impulse  responses  do  have  relatively  rapid  varia- 
tions (with  frequency)  compared  to  the  rates  of  variation  occurring  in  the  frequency  domain 
forms  of  most  test  (incident)  waveforms  commonly  used  in  TDR  or  TDT  measurements.  Conse- 
quently, the  application  of  the  homomorphic  deconvolution  transform  will  result  in  clustered 
groups  of  components  in  which  the  two  distinct  contributions  will  be  due  to  the  network 
impulse  response  and  the  incident  waveform,  respectively;  and  accordingly  they  will  be 
easily  separable.     Figure  3.2  shows  the  homomorphic  deconvolution  transformation  as  being 
applied  to  a  typical  time  domain  waveform  that  could  be  encountered  in  either  a  TDR  or  a 
TDT  measurement . 
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The  microwave  networks  that  will  be  considered  in  the  next  two  sections  of  this 
chapter  will  be  classified  into  the  following  two  classes : 

(1)  Pure  resistive  networks,  i.e.,  with  pure  resistive  reflections  (discontinuities), 
and  this  will  be  considered  in  section  3.3. 

(2)  Networks  with  reactive  discontinuities,  which  will  be  dealt  with  in  section  3.4. 

The  reason  for  this  classification  is,  as  will  be  seen  later,  that  a  pure  resistive  network 
has  a  train  of  delta  functions  as  its  impulse  response  in  comparison  with  reactive  reflec- 
tions, where  the  impulse  response  contains  a  train  of  pulses  whose  duration  and  time 
constants  increase  with  an  increase  in  the  reactive  parts  of  the  discontinuities.  This 
will  appear  in  the  cepstrum  domain  in  the  form  of  more  separability  for  the  resistive 
network  as  compared  to  reactive  ones. 

3.3    Pure  Resistive  Networks 

In  this  section  two  types  of  resistive  discontinuity  problems  will  be  considered: 
(1)  pure  resistive  loading  (discontinuities)  along  a  lossless  transmission  line,  figure  3.3 
and  (2)  the  case  of  commensurate  transmission  lines,  figure  3.4,  which  is  a  cascade  of 
lossless  transmission  lines,  each  one  having  different  electrical  parameters  (Z^,  y) .  In 
practice,  these  two  considered  cases  may  include  reactive  components  in  their  physical 
realization;  but  such  reactive  components  are  always  neglected  in  the  idealized  case,  and 
the  networks  are  treated  as  purely  resistive  ones.     However,  the  extension  to  cases  with 
more  general  reflections  including  reactive  parts  will  be  treated  in  the  next  section.  , 

In  this  section  the  time  domain  ref lectrometry  (TDR)  method  will  be  used.     The  general 
form  of  the  reflective  impulse  response  h'(t)  of  a  microwave  resistive  network  is  given  by 

h'(t)  =    I    A^6(t-T^)  (3.1) 
n=0 

where  t  is  the  time  variable  and  both  A    and  x     are  constants  determined  by  the  network 

n           n  ^ 

parameters.  Assuming  an  incident  waveform  x(t)  at  the  input  terminals  of  the  network,  the 
reflected  waveform  y'(t)  would  be  the  convolution  of  x(t)  and  h'(t);  i.e., 

y'(t)  =  x(t)  *  h'(t).  (3.2) 

The  observed  waveform  y(t)  at  the  input  gate  is  the  sum  of  both  incident  and  reflected 
ones ;  i.e., 

y(t)  =  x(t)  +  y'(t) 

=  x(t)  +  x(t)  *  h'(t) 

=  x(t)  *  (5(t)  +  h'(t)) 

=  x(t)  *  h(t)  (3.3) 
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where 

h(t)  =  6(t)  +  h' (t). 
In  the  frequency  domain  the  signal  h(t)  transforms  into  H(e-''^)  where 


HCeJ'^)  =  1  +    I  A 


e 
n=0 

■  00  '00 

=  1  +    5]    A    COS  0)  T    -  j         A    sin  w  x  (3.5) 
n=0    "  n=0    "  " 

which  is  clearly  seen  to  be  oscillatory  in  nature  (with  frequency).     On  the  other  hand 
X(e-''^)  for  most  physical  signals  is  non-oscillatory;  i.e.,  slowly  varying  with  frequency 
compared  with  H(e~''^).     The  application  of  the  homomorphic  deconvolution  can  then  be 
expected  to  separate  the  two  components  x(t)  and  h(t)  from  their  convolution. 

In  the  following  the  two  above  described  cases  of  resistive  networks  will  be  considered. 
In  each  case  two  examples  will  be  considered  and  developed;  the  first  of  which  will  be 
simple  (low  number  of  discontinuities)  such  that  the  mathematical  verification  of  the 
results  can  be  followed  without  too  much  confusion.     The  second  example  containing  a  large 
number  of  discontinuities  will  be  presented  in  terms  of  the  results  of  applying  the  homomorphic 
deconvolution  in  order  to  demonstrate  the  ability  of  the  method  to  handle  such  complicated 
cases. 

In  the  numerical  computations  for  the  examples  presented  below,   1024  samples  in  both 
time  and  frequency  domains  were  used.     The  time  and  frequency  sampling  intervals  were 
2.45  ps  and  0.4  GHz,  respectively. 

3.3.1    Case  of  pure  resistive  line  loading 

The  general  form  of  network  to  be  considered  in  this  part  is  given  in  figure  3.3. 

However,  as  mentioned  earlier,  only  a  simple  case  will  be  considered  as  far  as  the  mathe- 
matical details  are  concerned;  this  is  the  case  of  two  discontinuities  shown  in  figure  3.5. 

Two  resistances        and        are  connected  across  a  lossless  transmission  line  having  a 

characteristic  impedance  R^  at  the  distances        and  shown  in  the  figure.  Denoting 

the  reflection  coefficients  at  the  discontinuities  R^^  and  R2  as        and  respectively, 

the  reflective  impulse  response  of  the  network  is  derived  in  figure  3.5b  to  be 

h'(t)  =  p^6(t-2T^)  +  p^(l+p^)^  6(t  -  (21^+2X2))  +  ....  (3.6) 

Upon  substituting  for  p^,  P2,  x^,  and  X2  by  their  respective  values  -0.2,  -0.75,  127  ps, 
and  84.67  ps  that  were  used  in  the  simulated  application,  this  impulse  response  becomes 

h'(t)  =  -0.2  6(t-3x)  -  0.48  6(t-5x)  -0.072  6(t-7x)  -    (3.7) 
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where  t  =  84.67  ps .     Figure  3.5c  shows  the  impulse  response  h(t)  =  6(t)  +  h' (t)  for  this 
specific  case. 

Figure  3.6  shows  the  time  domain  waveform  Sg(t)  of  the  generator  used,  while  figure  3.7 
shows  the  corresponding  frequency  domain  form.     The  incident  waveform  x(t)  is  given  by 
e  (t) 

x(t)  =  -^2~  (3.8) 

since  Rq  the  generator  internal  resistance  is  matched  to  the  line  resistance.    Thus  x(t) 
can  also  be  represented  by  the  figures  3.6  and  3.7  scaled  by  a  factor  of  one  half. 

Figures  3.8  through  3.14  show  various  processing  waveforms  throughout  the  steps  of 
the  homomorphic  deconvolution  transform.     Figures  3.8  and  3.9  show  the  observed  waveform 
at  the  input  gate  of  the  network  in  both  time  and  frequency  domains,  respectively.  Figure 
3.10  shows  the  resultant  cepstrum  and  the  linear  filters  to  be  used  for  the  separation 
process.     Figures  3.11,  3.12,  3.13,  and  3.14  give  the  time  and  frequency  forms  of  the 
deconvoluted  components  x(t)  and  h(t),  respectively.     The  obtained  forms  of  x(t)  and  h(t) 
do  coincide  with  the  original  form  to  a  great  extent.     The  erroneous  oscillations  observed 
in  the  deconvoluted  impulse  response  is  in  fact  due  to  sampling  errors;  i.e.,  finite 
number  of  samples  and  finite  observation  time  (window) . 

The  figures  3.15  through  3.19  show  another  example,  but  a  more  complicated  one  with 
six  discontinuities  along  the  line.     Figure  3.15  shows  the  network  considered  and  various 
values  of  various  network  parameters.     The  generator  signal  used  is  the  same  one  as  given 
in  figures  3.6  and  3.7.     Figure  3.16  shows  the  convolution  y(t),  while  figure  3.17  shows 
the  corresponding  cepstrum.     The  deconvoluted  x(t)  and  h(t)  are  given  in  figures  3.18  and 
3.19,  respectively. 

3.3.2    Commensurate  transmission  lines 

The  general  form  of  a  commensurate  transmission  line  network  was  given  in  figure  3.4. 
Again,  two  examples  will  be  considered  in  this  part,  one  example  with  two  discontinuities 
and  the  other  one  with  six  discontinuities.     The  same  generator  signal  that  was  used  in 
section  3.3.1  will  be  used  in  this  section. 

The  two  discontinuity  example  is  shown  in  figure  3.20.     Figure  3.20a  gives  the 
network;  it  is  a  cascade  of  two  lines  having  the  following  parameters: 
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Characteristic 
Impedance 


Line  1  Line  2 

\  +  jO  +  jO 


2 


Propagation  n  ■    •    /  r.  ■    .  / 

Constant  -^1 

Length     .  ^1  ^2 

s  One  Way  „    ,  „  , 

•  Travel  Time  ^1  =  V^l  ^2  =  V^2 

The  exiting  generator  has  a  matching  impedance  R^.     Figure  3.20b  shows  the  multiple 
reflections  pattern.     The  reflective  impulse  response  as  can  be  written  from  this  figure  is 

h'(t)  =  p6(t-2T^)  -  p(l-p2)  6(t  -  (2x^+2x2))  -  p3(l-p2)  6  (t  -  (2x^+4x2))  -  ....  (3.9) 

While  h(t)  =  6(t)  +  h' (t)  is  shown  in  figure  3.20c  for  the  specific  case  that  is  to  be 
simulated  on  the  computer,  in  which 

p  =  0.333,  .  , 

x^  =  127  ps, 

and 

X2  =  101.6  ps,        ■     ..  '  ' 

this  h(t)  is  given  by 

h(t)  =  6(t)  +  0.333  6(t-3x)  -  0.2963  6(t-5.4x)  -  0.0329  5(t-7.8x)  -    (3.10) 

Figures  3.21  through  3.27  show  the  processing  waveforms    throughout  the  steps  of  the 
homomorphic  deconvolution  transform.     Again,  figures  3.25  and  3.27  give  the  deconvoluted 
x(t)  and  h(t)  which  are  demonstrated  to  be  in  complete  agreement  with  the  original  forms 
for  x(t)  and  h(t).     On  the  other  hand  figures  3.28  through  3.32  show  the  six  commensurate 
lines  case. 

3.4    Networks  With  Reactive  Components 

In  the  present  section  two  additional  applications  will  be  introduced  beyond  the 
generalization  to  more  general  discontinuities  along  the  line.     Section  3.4.1  will  be 
devoted  to  the  consideration  of  a  TDT  type  measurement,  while  section  3.4.2  will  deal  with 
the  effects  of  an  unmatched  source  impedance  in  a  conventional  TDR  measurement. 
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3.4.1    Time  domain  transmission  (TDT) 

Figure  3.1b,  earlier  in  the  chapter,  showed  the  experimental  setup  for  a  typical  TDT 
measurement.     The  example  that  will  be  considered  in  this  section  is  the  following: 
two  impedance  discontinuities        and        were  connected  across  a  line  having  characteristic 
impedance        and  a  propagation  constant        at  distances  d^  and        as  shown  in  figure  3.33. 
A  generator  of  matched  impedance  excites  the  network;  ^g(t)  is  the  generator  time  domain 
waveform.     The  aim  of  this  specific  problem  was  to  see  how  the  TDT  waveforms  would  yield  to 
the  homomorphic  deconvolution  transform. 

In  the  numerical  computations  for  the  present  example,   1024  samples  in  both  time  and 
frequency  domains  were  used.     The  time  and  frequency  sampling  intervals  were  4.9  ps  and 
0.2  GHz,  respectively. 

The  values  of  the  network  parameters  that  were  used  in  this  specific  example  were  as 
follows : 

line:     7  mm  precision  coaxial  line  (50  H)  [6] 
Z^:     Capacitive,        =  0.3  pf 
Z^:     Capacitive,        =  0.1  pf 
d2 :     3.54  cm 

The  waveform  of  the  generator  ^gCt)  is  given  in  figure  3.34.     The  network  was  simulated  on 

a  digital  computer,  and  the  resultant  TDT  waveform  at  the  outout  terminals  e  (t)  is  given 

o 

in  figure  3.35.     The  frequency  domain  form  E^Ce-'"')  of  e^(t)  is  given  by  the  relation 

E  (e^^)  exp[-Y^(o:)- (d  +d  )]  [1+p   (O]  [l+c,(oj)] 
°  2  [l-p^(ai)P2(!^)expr-2Y^(uj)d2)] 

where  E  (e-^'^)  is  the  frequency  domain  form  of  e  (t),  and  c^C^)  and  c-C-)  are  the  reflection 

g  g  1 

coefficients  corresponding  to  the  impedance  discontinuities  Z^(^)  and  Z^{^) ,  respectively. 

The  magnitude  of  E^Ce-'")  is  given  in  figure  3.35.     The  homomorphic  deconvolution  transform, 

figure  3.2,  is  to  be  applied  on  the  output  signal  to  separate  its  convolution  components. 

As  was  mentioned  before,  the  filtration  in  the  cepstrum  domain  separates  two  components 

that  have  substantially  different  rates  of  variation  in  their  frequency  domain  forms.  By 

examining  eq,    (3.11)  it  can  be  concluded  that  the  slowly  varying  component  E     (e"'^)  would 

OS 

be 

E  (ej") 

E     (e^")  =    exp[-7  (^)-(d^+d„)]     [1+P,(-.)]   [l+=,(^)]  (3.12) 

OS  ^  o  L     I  1  I 

which  is  the  first  path  signal  through  the  network,  while  the  rapidly  varying  one  ^^^(e-'") 
is  the  multiple  reflection  effect, 

E     (e^'^)  ■   .  (3.13) 

or 


[l-p^((jo)p2(w)  ]  exp  -2Y^(aj)d2 
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These  eqs.    (3.12)  and  (3.13)  can  be  used  to  determine  up  to  two  unknown  system  parameters, 
provided  that  all  other  network  parameters  are  known.     To  show  this,  eq.    (3.13)  will  be 
used  to  calculate  P2(<Jj)  assuming  that  p^(e-''^)  is  known;   the  values  of  P2((jj)  are  to  be  used 
to  calculate  the  corresponding  impedance  Z2(w). 

The  above  mentioned  steps,  in  the  current  section,  were  followed,  and  the  results  are 
shown  in  figures  3.37  through  3.42.     Figure  3.37  gives  the  cepstrum  of  e^(r).     Figures  3.38 
through  3.41  give  the  slowly  and  rapidly  varying  components  in  both  frequency  and  time 
domain  forms,  respectively.     The  resultant  values  for  E^^(e-'^)  were  used  for  the  calculation 
of  Z^iiii)  as  indicated  earlier.     The  computed  values  for  Z^(ii>)  showed  a  pure  capacitive 
impedance,  and  the  capacitance        is  plotted  in  figure  3.42.     The  figure  shows  erroneous 
oscillation  about  the  actual  given  value.     These  error  oscillations  are  again  due  to 
sampling  errors.     Part  of  the  errors  are  also  due  to  the  imperfect  separation  in  the 
cepstrum  domain  because  of  the  extended  long  tails  due  to  reactive  reflections  on  the  line. 

3.4.2     Time  domain  reflections  (TDR)  with  reactive  loading/unmatched  generator  impedance  [7,8] 

In  this  case  the  TDR  measuring  technique  will  be  presented  in  its  conventional  form, 
showing  the  limitations  of  the  method.     Then  it  will  be  shown  how  the  homomorphic 
deconvolution  transform  helps  in  offsetting  these  limitations. 

3.4.2.1     TDR  measurements 


Figure  3.43  shows  the  typical  TDR  experimental  setup  for  measuring  the  reflection 
coefficient  of  the  line  termination  p^(aj).  The  general  procedure  is  to  record  the  dis- 
played signal  for  the  termination  under  test  and  two  other  reference  terminations,  usually 
a  short  circuit  and  a  Z^(a))  termination.*  The  complex  frequency  domain  forms  of  the  three 
recorded  signals  are  given  by 


Z  (o)) 

E^^(eJ'')  =  E  (e^")  •   ^   •  T((.)  (3.14) 


dM' 


Z  (co)  +  Z  (uj) 
o  e 


E,.(ej'^)  =  E  (ej'^)  1  -  exp  f-2Y  (o))   ^^^^  (3_^5) 


Z  (to)  +  Z  (w)       1  +  p   (w)  •exp(-2Y(a3)£} 
o  g  g 

Z(i^)                 1  +  p   (to)  exp(-2Y(u))£} 
E^Ae^"")  =  E  (eJ'^)  ^  h  T(to)  (3.16) 


dL' 


Z  (03)  +  Z  (o))  1  -  P,((o)p  (to)  exp  f-2Y(to)£) 
o  e  L        "  '  ^ 


corresponding  to  Z  (to),  short  circuit,  and  Z  (tjj)  terminations,  respectively. 

O  Li 

One  cannot  solve  the  above  three  equations,  as  they  are,  to  obtain  p^((d), 


*For  an  exciting  waveform  which  becomes  of  absolutely  constant  value  before  the  reflections 
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arrive,  the  Z    measurement  is  not  necessary, 
o 


3.4.2.2     Conventional  nethod 

The  conventional  method  for  evaluating       (-)  is  to  consider  only  the  TDR  signals 
occurring  in  the  tirae  window  0  <  t  <  41,  where  I  is  the  travel  time  on  the  line  of  length 
I.     In  this  case  ecs.    (3.11),    (3.12),  and  (3.13)  will  take  the  simplified  foms 


"dM,4T 


Z  (_) 

(eJ*^)  =  E  CeJ'")   •   5   .  T(cj)  (3.17) 

^  Z^(_)  +  Z(^) 

o  g 


Z  (:j) 

E^c  ATfe-"")  =  S  (e^'')  •   •  T(-)  •   [1  -  (l-H)„(«a))  exp (-2y (o) l)  ]  (3-18) 

"^'^^  S  Z  (O  -  2   (.)  S 


7  C 


E..  ^^ie^"")  =  E^(e^-)   ■   ^   •  T(u)   •   [1  +  p.  (to)   [l-^pC^)]  exp  (-2y  (w)     ]  (3.19) 

""^'"^  ^  Z  (-)  ^  Z   (-)  ^  ^ 

o  g 

fron  which       (_)  can  be  deteriiined; 

E,.   ,_(-■)  -  E^,^,  _(w) 

pjcj)  =  -   ^=^^^^^^         .  (3.20) 

This  tiethod  vorks  only  if  the  second  and  subsequent  reflections  do  not  overlap  the 
first  one.     Figure  3.44  shows  two  cases  of  TDR  signals:     the  non— overlapping  reflezticns 
where  the  conventional  method  works,  and  the  overlapping  reflections  in  which  the  conventional 
method  fails . 

3.4.2.3    The  application  of  the  ho-oncrphic  ceconvolution  method 

In  the  homomorphic  deconvolution  method,   see  figure  3.2,   the  ccnplete  forr.s,  ecs. 

(3.11),   (3.12),  and  (3.13),  are  used  to  obtain  not  only       (^)  but  also  c^(!j)),  7(w),  and  the 

product  E^(e^'^)T(ii))  .     The  process  is  begun  by  first  obtaining  the  differences  of 
eq.    (3.15)  -  eq.    (3.14)  and  eq.    (3.16)  -  eq.    (3.14)  which  are 

^S^^'")  =  ^ds(^''>  -  ^dM^^''^ 

Z  (t^) 

=  E  (e^")  •   •  T(.)   •  exp{-2v(M)£) 

S  Z  ('^)  -  Z  (_) 

O  3. 


(-1)    •     1  ^  


:1  -r  z^U)  eKp[-2y(m)i] 


(3.21) 
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=  E  (e-^^) 


Z(i^)  +  Z  (o)) 
o  g 


T(u))   •  exp(-2Y(a))JlJ 


1  +  P_(w) 


1  -  p   (w)p^(a))  exp[-2Y(a))£)J 


Using  the  homomorphic  deconvolution  technique,  each  of  the  two  equations,  eqs.  (3.21) 
and  (3.22),  can  be  split  into  two  components,  the  slowly  varying  component  and  the  rapidly 
varying  one  (subscripted  by  s  and  r  in  the  following  equations,  respectively). 


E_  (eJ-^)  =  E  (eJ") 
Ss  g 


Z  (cj)  +  Z  (u) 
o  g 


T(a))   •  exp(-2Y(tj)£)    •  (-1) 


1  +  P  (oo) 


1  +  P   (w)-exp(-2Y  (w)«.] 


(e-J  )  =  E  (e-J  ) 
Ls  B 


Z  (w)  +  Z  (w) 

O  2 


T(a))   •  exp(-2Y(a)))l)    •  p^(a)) 


1  +  P_(a)) 


1  -  p  (w)'Pj^(u))'exp(-2Y(a))X.) 


The  above  four  equations  together  with  eq.    (3.14)  can  be  used  to  solve  for  p   (w),  p   (oi) , 
iw  L  g 

y(oj),  and  the  product  Eg(e     )T(a)) . 

The  homomorphic  deconvolution  technique  applied  to  TDR  problems  not  only  gives  access 

to  more  of  the  microwave  network  parameters  but  also  works  for  overlapping  reflection 

cases,  where  the  conventional  methods  fails. 


3.4.2.4  Application 


A  computer  program  was  written  to  simulate  the  TDR  network  of  figure  3.43b  with  parallel 
R-C  combinations  to  serve  as  reflective  load  and  generator  impedances.     The  homomorphic 
deconvolution  technique  was  applied  to  the  waveforms  of  the  simulated  model;  the  network 
parameters  and  various  waveforms  are  given  in  figures  3.45  through  3.50. 

In  the  numerical  computations  for  the  present  example  1024  samples  were  used.  The 
time  and  frequency  sampling  intervals  were  2.45  ps  and  0.4  GHz,  respectively. 

Equations  (3.23)  and  (3.25)  were  used  to  calculate  the  reflection  coefficient  Pj^(w) 
from  the  transform  domain  data  of  the  separated  signals  (shown  in  the  time  domain,  figures  3.48 
and  3.50),     Then  the  resultant  Pj^(<^)  was  used  to  determine  R  and  C  of  the  load  impedance 
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Z  (o)) .     The  calculated  values  of  R  and  C  vs.  frequency,  figures  3.51  and  3.52,  oscillate 

Li 

around  their  respective  actual  values.     These  erroneous  oscillations  which  correspond  to 
those  indicated  in  figures  3.48  and  3.50  are  again  due  to  (a)  the  finite  number  of  samples 
considered,   (b)  the  finite  observation  time  window,  and  (c)  the  imperfect  separation  in  the 
cepstrum  due  to  overlapping  long  tails  of  the  different  components. 

3.5    Summary  and  Conclusions 

The  present  chapter  presented  several  classes  of  time  domain  measurement  problems  in 
the  field  of  linear  microwave  network  analysis.     The  homomorphic  deconvolution  transform, 
that  was  developed  in  chapter  2,  was  used  in  deconvo luting  various  convolutions  in  these 
problems.     The  examples  that  were  considered  are  as  follows: 

(1)  Pure  resistive  loading  along  a  transmission  line. 

(2)  Commensurate  transmission  line  networks. 

(3)  TDT  type  measurement  with  two  reactive  discontinuities. 

(4)  Conventional  TDR  measurement  to  determine  a  reactive  load  impedance. 

The  homomorphic  deconvolution  method  was  proven  to  be  very  effective  in  separating 
convolutions.     The  method  was  not  limited  by  the  overlap  of  consequent  reflections  in  the 
time  domain. 
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a.  NETWORK 


eg(t)l 


TIME 
0 


2t, 


SIGNAL  LEVEL 
1 


2  2 
2t^+4t2     P1P2  (l+P-]) 

2t^+6t2     P^^P2"^(1+P]  )^ 

2t^+8t2  P^"^P2^(1+P] 
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-.1 
.  -.2 
-.3 
-.4 
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h(t)  =  6(t)+h'(t) 


.2,  p2 


=  -.75 


1   '  '2 

WHERE  I  =84.67  ps 
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Figure  3.5.    Case  of  Two  Reflections. 
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Figure  3.20.   Case  of  Two  Commensurate  Transmission 
Lines  with  Two  Discontinuities. 
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Figure  3.33.    TDT  Network 
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Figure  3.50.         (t) ,  Arrows  Indicate  Error  Oscillations. 
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4.     TIME  DOMAIN  DIELECTRIC  SPECTROSCOPY 


4.1  Introduction 


In  chapter  3  the  homomorphic  deconvolution  transform  was  successful  in  solving  various 
time  domain  network  analysis  problems,  and  more  specifically,  time  domain  ref lectometry  and 
transmission  problems.     Network  parameters  of  interest  were  identified,  and  more  importantly, 
the  time  domain  overlap  of  multiple  reflections  did  not  put  any  limitations  or  restrictions 
on  the  method. 

In  the  present  chapter  another  application  of  time  domain  measurements  will  be  con- 
sidered.    Here,  the  homomorphic  deconvolution  transform  will  be  applied  to  time  domain 
dielectric  spectroscopy  (TDDS);  i.e.,  the  study  of  dielectric  relaxation  phenomenon  by  time 
domain  methods. 

TDDS  is  in  fact  a  TDR  problem,  but  with  the  emphasis  being  turned  to  the  dielectric 
material  filling  the  transmission  line  in  order  to  evaluate  the  frequency  dependent  complex 
dielectric  constant  (relative  permittivity)  k(w) . 

As  was  indicated  in  chapter  3,  the  conventional  TDR  method  is  based  upon  and  limited 
by  time  windowing  to  select  the  relevant  or  desired  waveform  epoch.     Thus,  for  the  con- 
ventional TDR  method  to  be  used,  the  dielectric  sample  has  to  be  either  relatively  infini- 
tesimal or  relatively  infinite  in  length  to  obtain  nonover lapping  time  windows.  The 
homomorphic  deconvolution  method  will  be  applied  in  the  present  chapter  to  remove  these 
restrictions  and  allow  the  use  of  finite  size  samples. 

The  method  will  be  applied  to  a  model  typical  of  dielectric  materials  having  polariza- 
tion absorptions  over  the  frequency  range  from  1  kHz  upward  to  100  GHz. 


The  typical  experimental  setup  for  a  TDDS  measurement   [12-18]  is  shown  in  figure  4.1. 
The  dielectric  material  to  be  characterized  is  formed  into  the  shape  required  to  fill  a 
coaxial  sample  holder  of  finite  length  d  corresponding  to  an  empty  precision  line  having  a 
characteristic  impedance  Z^(w)  and  a  propagation  constant  y(<^)'     A  time  domain  waveform 
e^(t)  is  to  be  recorded.     This  waveform  consists  of  the  multiple  reflections  from  both 
dielectric-air  interfaces  as  well  as  the  incident  wave. 

Because  the  physical  dimensions  of  the  two  transmission  lines  are  the  same,  the  charac- 
teristic impedance  Z^(a))  and  the  complex  propagation  constant  Yj(w)  of  the  dielectric 
filled  transmission  line  are  related  to  the  corresponding  quantities  of  the  uncelled  air 
line  as  follows: 


4.2    Theory  of  Time  Domain  Dielectric  Spectroscopy 


Z  (w)/A:(w) 


o 


(4.1) 


and 


Y  (u)    •  vk(oo) 


o 


(4.2) 
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where  k(a))  is  the  complex  dielectric  constant  of  the  sample.  The  reflection  coefficient  at 
the  air-dielectric  interface  on  the  line  p (w)  is  given  by 


p(oo)  =  [Z^(a)) 


-  Z  (w)]/[Z,(ai)  +  Z  (o))] 
o  a  o 


(4.3) 


which  upon  substituting  eq.    (4.1)  takes  the  form 


pM  =  [l-/k(w)]/[H-/k(a))]. 


(4.4) 


The  form  of  eq.    (4.4)  can  be  solved  for  k(oj)  as  follows: 


A(a3)  =  (l-p(a)))/(l+p(cj)) 


br 


k(a))  =  [(l-p(a)))/[l+p(oo))]2. 


(4.5) 


This  equation,  eq.    (4.5),  summarizes  the  concept  of  using  the  TDR  method  in  dielectric 
characteristics  by  experimentally  observing  the  reflected  wave  produced  by  the  reflection 
coefficient  p (w) .     The  deduced  reflection  coefficient  p(a))  is  then  substituted  in  eq.  (4.5) 
to  yield  the  complex  dielectric  constant  k(a)). 

In  the  conventional  TDDS  method  [12-18]  the  first  reflection  coming  back  from  the  air- 
dielectric  interface  determines  the  measurement's  time  window;   i.e.,  the  reflected  waveform 
is  observed  for  the  time  interval  containing  the  first  reflection.     This  first  reflection 
has  a  Fourier  transform  which  is  the  product  of  the  Fourier  transforms  of  the  incident 
signal  e  (t)/2  and  the  reflection  coefficient.     The  separability  of  the  time  window  and  the 
knowledge  of  ^g(t)  allow  the  calibration  of  p(a))  and  ultimately  k(u) . 

As  mentioned  in  the  introduction  section,  the  separability  of  the  first  reflection 
requires  either  an  infinitesimal  or  an  infinite  sample  length.     The  requirement  of  having  a 
very  long  sample  is  not  generally  practical.     On  the  other  hand  an  infinitesimal  sample 
would  not  alter  ^g('^)  ^°  great  degree,  thereby  not  imparting  to  the  signal  waveform 

much  information  about  the  dielectric.     A  finite  sample  of  practical  size  will  produce 
overlapping  multiple  reflections  thus  invalidating  the  conventional  (windowing)  method.  In 
such  cases,  the  homomorphic  deconvolution  transform  will  succeed  in  removing  (deconvolu- 
ting)  the  effects  of  multiple  reflections  allowing  the  determination  of  p(co). 

The  frequency  domain  form  E  (e-''^)  of  the  recorded  waveform  e  (t)  is  given  by 


By  virtue  of  eqs.    (4.2)  and  (4.4),  the  above  eq.    (4.6)  has  only  one  unknown  k(a))  provided 

that  E  (e-''^),   Z  (w),  and  y  M  are  known.     However,  such  an  equation  can  not  be  solved,  as 
go  o 

it  is,   to  yield  this  unknown.     Some  research  workers   [12]  have  gone  through  approximations 
to  simplify  this  equation;  but  the  use  of  the  homomorphic  deconvolution  transform  on 
eq.    (4.6),  as  it  is,  will  yield  the  required  information  needless  of  any  assumptions  or 
approximations. 


1  -  exp  (-2Yj (w)d) 


(4.6) 
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4.3    The  Homomorphic  Deconvolution  as  Applied  to  TDDS  Problems 

The  homomorphic  deconvolution  transform  that  was  first  developed  in  chapter  2,  and 
then  used  with  TDR  networks  in  chapter  3,  separates  a  convolution  into  its  two  components. 
These  two  components  are  the  slowly  and  rapidly  varying  frequency  domain  components.  By 
examining  eq.   (4.6)  the  slowly  and  rapidly  varying  components  can  be  identified  as 

E  (e>) 

(e^")  =-^-5         exv[-2y(ui)l]   p  (oo)   ,  (4.7) 

Is  /  o 

and 

1  -  exp  f-2Y  ,  (co)dl 

=  ■   (4.8) 

^  1  -  p^(a))  exp  (-2Y^(tD)d) 

respectively.     In  these  two  equations,  eq.   (4.7)  is  a  simple  one  that  can  be  solved  in 

p(u)),  and  consequently  k.(to) ,  while  eq.   (4.8)  still  has  the  same  complications  that  originally 

existed  in  eq.  (4.6). 

4.4    Model  for  Debye  Dielectric  [9,  11] 

One  well  known  class  of  dielectric  materials  will  be  considered  in  the  application 
presented  in  this  chapter,  the  single  relaxation  time  dielectric  known  as  the  Debye  dielec- 
tric.    The  complex  dielectric  constant  k(a))  of  a  Debye  dielectric  is  given  by 

^    ~  K 

k(w)  =        +  —   (4.9) 

1   +  jtOT 

where  k    =  k(to=0) ,  k    =  k(u='»)  ,  and  t  is  the  dielectric  relaxation  time.     This  equation  is 

O  CO 

commonly  put  in  the  form 

k(w)  =  k'(uj)  -  jk"(ca)  (4.10) 

where 

k'(aj)   =    (k^+uj2^2i^^-)/d+^2^2)  (4.11) 

and 

k"(u)  =  wT(k  -k  )/(l-Ko2^2)_  (4.12) 

O  CO 

The  plot  k"(w)  vs.  k'(u)  is  a  semicircle  of  radius  r  and  centered  at  the  point  (c,0) 

where 

r  =  (k^-kJ/2  (4.13) 

and 

^  =  (V^-)/2.  (4.14) 
the  point  (c,r)  is  the  point  at  which  to  =  u^,  where 

V  =  -81-  (^-15) 


The  semicircle  is  shown  in  figure  4.2;  this  semicircle  is  commonly  called  the  Debye  semi- 
circle or  the  Cole-Cole  plot  [10-11]. 

It  is  important  to  mention  that  the  complex  dielectric  constant  parameters  k^,  k^,  and  t 
can  be  determined  once  the  dielectric  semicircle  is  obtained.     This  is  achieved  by  virtue 
of  the  eqs.    (4.13),   (4.14),  and  (4.15)  to  be 

k    =  c  +  r, 
o 

k    =  c  -  r, 

00 

and 

T  =  l/oj  . 
o 

Also  the  circular  geometry  of  the  graph  provides  a  means  of  consistently  smoothing  experi- 
mental data  for  k^(a))  and  k'"(aj). 

Finally,  the  point  should  be  made  that  the  TDDS  methods  developed  here  and  not  restricted 
to  any  particular  dielectric  model.     The  Debye  model  has  been  chosen  because  it  is  a  well 
known  model  for  the  dilute  mixture  of  polar  molecules  in  nonpolar  solvents.     As  such  the 
model  will  clearly  demonstrate  the  utility  of  the  deconvolution  method  in  obtaining  the 
dielectric  parameters,  k''(aj)  and  k"(aj). 

In  general,  a  dielectric  material  may  contain  a  distribution  of  relaxation  times 
comprising  its  relaxation  process.     The  relative  dielectric  constant  may  be  viewed  as  a  sum 
of  relaxation  processes  each  having  a  distribution  of  relaxation  times  [19]. 

oo  , 

k(a))  =  k    +  /    dT 

o 

where  g(T)  is  the  relaxation  time  distribution  function. 


4.5    Application  of  Homomorphic  Deconvolution  in  TDDS  [20] 

A  computer  program  was  written  to  simulate  the  TDDS  setup  of  figure  4.1.     The  transmis- 
sion line  characteristic  impedance  was  50  U  and  the  dielectric  sample  length  was  one  cm. 
The  dielectric  parameters  were  chosen  to  by  typical  of  actual  dielectrics  satisfying  the 
Debye  model.     The  relative  dielectric  constant  employed  was 

^.r  ^      /   «  J.      24.5  -  4.8 
k(a))  =  4.8  +   


1  +  j2.0xl0  ^^w 


in  which  k    =24.5,  k    =4.8,  and  t  =  2  ps.     The  time  domain  waveform  of  the  generator 

o  =° 

e  (t)  is  the  same  one  used  in  section  3.4.2  of  chapter  3;  it  is  given  in  figure  3.46. 
S  -1 

The  homomorphic  deconvolution  transform  was  applied  to  the  computed  signal  E^(e-^   )  to 

separate  both  its  slowly  and  rapidly  varying  components.     The  slowly  varying  component 
together  with  eq.    (4.7)  were  used  to  calculate  the  reflection  coefficient  p(a)).     The  result- 
ant values  of  p(a))  were  substituted  in  eq.    (4.5)  to  yield  the  desired  k(uj) .     Figures  4.4, 
4.5,  and  4.6  show  the  results  of  the  above  mentioned  steps.     Figure  4.4  gives  the  magnitude 
of  E^(e-'*^);  figure  4.5  gives  the  magnitude  of  the  separated  slowly  varying  component  while 
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figure  4.6  shows  the  Cole-Cole  plot  of  the  calculated  complex  dielectric  constant  components 
(Argand  diagram).     A  semicircle  was  fitted  [21]  to  the  computed  plot  and  also  is  given  in 
the  figure.     The  cemicircle  parameters  were  r  =  9.80255,  c  =  14.60045,  and        =  0.50071  x 
10''"^  rad/sec. 

Upon  substituting  these  values  in  eqs.   (4.16),   (4.17),  and  (4.18),  the  dielectric 
parameters  are  evaluated  to  be 

k  =24.403, 
o 

k    =  4.7979, 

CO 

and 

T  =  1.997  ps, 

which  are  within  0.3%  error  of  their  originally  given  values. 

4.6    Summary  and  Conclusions 

In  this  chapter,  the  theory  of  the  time  domain  dielectric  spectroscopy  measurements 
was  presented.     Conventionally  this  method  of  measurement  works  only  for  non-overlapping 
reflections;  i.e.,  reflections  occurring  in  separate  time  domain  windows.     If  the  reflec- 
tions overlap,  only  approximate  solutions  may  be  applied  to  the  frequency  domain  form  of 
the  signal. 

By  applying  the  homomorphic  deconvolution  transformation  to  the  TDDS  waveform,  a  new 
method  was  developed  which  was  not  limited  by  overlapping  reflections.     The  new  method  was 
applied  to  the  well  known  Debye  dielectric  model.     The  method  eliminated  the  overlapping 
reflections  and  determined  the  dielectric  parameters  to  within  0.3%  error  of  their  model 
values . 

The  new  method  is  a  general  one  and  has  a  great  utility  for  time  domain  measurements 
on  dielectrics  since  it  does  not  include  any  restrictions  on  the  dielectric  material 
parameters  or  require  the  assumption  of  a  prior  physical  model.     Furthermore,  the  new 
method  is  not  restricted  to  time  domain  spectroscopy  of  only  dielectric  materials  but  has 
utility  for  time  domain  spectroscopy  in  general  such  as  for  magnetic  and  other  phenomena. 
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5.     THE  APPLICATION  OF  THE  HOMOMORPHIC  DECONVOLUTION  TRANSFORM  TO  THE 
REMOVAL  OF  SCATTERING  COMPONENTS  FROM  ANTENNA  TEST  RANGE  DATA 

5.1  Introduction 

Time  domain  analysis  of  electromagnetic  radiation  problems  has  become  increasingly 
more  popular  in  the  recent  years  because  the  time  domain  method  provides  physical  insight 
into  the  behaviour  of  antenna  structures  excited  by  pulsed  sources.     Time  domain  measure- 
ments do  offer  significant  advantages  over  the  classical  steady-state  frequency  domain 
measurement;  they  are  as  follows: 

a.  There  is  no  need  for  an  anechoic  chamber  to  absorb  unwanted  reflections  because 
time-windowing    can  be  shut-out  reflections. 

b.  There  is  no  need  for  repeated  magnitude  and  phase  measurements  at  different  fre- 
quencies because  the  time  domain  waveform  contains  the  complex  spectral  information  in 
a  single  entity. 

Frequency  domain  information  can  be  obtained  from  the  measured  time  domain  waveform  by 

simply  computing  the  corresponding  Fourier  transform.     However,  a  recently  applied  method 

[22,  23  and  24]  based  upon  Prony's  algorithm  [25]  can  also  be  used  to  deduce  the  required 

frequency  domain  information.     The  method  is  based  on  the  fact  that  the  time-dependent 

function  can  be  expressed  as  the  sum  of  complex  exponentials 

a.t  • 
f(t)  =  ^  A.  e  ^  (5.1) 

i 

for  which  the  frequency  domain  form  is 
A. 

v-^      '  ^    ■         ■ .  \      ;  , 

where       are  the  complex  poles  of  F(s)  and  A^  are  the  corresponding  complex  amplitudes  or 
residues.      By  fitting  the  time  domain  data  to  the  form  of  eq .   (5.1)  and  determining  the 
location  of  the  singularities  a^'s  and  the  corresponding  residues  A^'s,  the  system  under 
measurement  is  fully  characterized  according  to  classical  circuit  theory. 

This  method  has  the  advantage  of  yielding  an  analytical  spectral  response  for  the 
experimental  system  rather  than  the  discrete  form  that  the  numerically  implemented  Fourier 
transform  yields.     Also,  it  is  capable  of  extrapolating  the  measured  transient  response  for 
both  early  and  late  times. 

It  is  clear  that  the  above  discussion  favors  time  domain  measurements  for  pulsed 
electromagnetic  radiating  systems.     Also,  it  is  clear  that  it  should  be  possible  to  window- 
out  any  unwanted  reflections  and  scattering  components.     However,  for  cases  where  these 
undesired  reflections  or  scattering  components  do  overlap  on  the  useful  part  of  the  time 
domain  signal,  the  homomorphic  deconvolution  should  enable  the  separation  and  the  rejection 
of  the  unwanted  components. 
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The  present  chapter  will  demonstrate  the  power  of  the  homomorphic  deconvolution  trans- 
form in  separating  and  removing  unwanted  reflections  and  scattering  components  in  a  pulsed 
electromagnetic  radiating  system.     The  transform  will  be  applied  to  radiation  from  a  wide- 
band horn  antenna  pulsed  by  baseband  signal  having  frequency  components  out  to  beyond 
6  GHz. 

5.2    The  Homomorphic  Deconvolution  as  Applied  to  Antenna  System 

The  homomorphic  deconvolution  transform  was  developed  earlier  in  chapter  2,  and  was 
used  in  chapters  3  and  4  for  the  separation  of  multiple  reflections  due  to  different  dis- 
continuities on  a  transmission  line.     In  these  cases  point  or  discrete  discontinuities 
generated  reflections  which  transformed  to  clusters  (discrete  well-separated  pulses)  in  the 
cepstrum-domain.     The  localized  clusters  were  then  separated  by  filtering. 

In  radiation  systems,  the  situation  is  different.     In  addition  to  the  point  reflections 
there  exist  plane  reflections,  which  can  be  called  a  continuium  of  multipaths.     Also,  there 
exist  scattering  components  due  to  the  interfaces  between  the  antenna  aperture  and  the 
radiation  medium,  and  also  those  due  to  different  nearby  scattering  objects.     Figure  5.1 
illustrates  some  of  these  phenomena. 

Consequently,  it  is  not  expected  that  definite  or  point  located  clusters  will  appear 
in  the  cepstrum-domain  of  such  systems;  rather,  it  is  expected  that  a  distribution  of 
contributions  will  appear  in  the  cepstrum-domain.     Figure  5.2  shows  a  typical  cepstrum  for 
antenna  systems.*    This,  of  course,  makes  the  separation  of  unwanted  components  more  dif- 
ficult than  the  cases  of  chapters  3  and  4.     Accordingly,  the  degree  of  filtering  must  be 
determined  by  a  trial  and  error  procedure.     If  too  much  of  the  cepstrum  spread  components 
is  removed,  it  may  include  parts  of  the  useful  signal  information.     On  the  other  hand,  not 
removing  enough  of  the  unwanted  components  may  result  in  a  remainder  which  overlaps  the 
useful  time  domain  signal. 

The  separation  in  the  cesptrum  domain  then  should  be  done  on  a  successive  trial  basis 
using  different  filters;  the  results  should  be  examined  in  the  light  of  the  physical  system 
to  decide  upon  the  best  filtering  process. 

The  next  section  of  the  present  chapter  will  be  devoted  to  the  description  of  the 
antenna  system  to  be  used  in  the  forthcoming  application. 

5.3    The  Antenna  System 

The  experimental  setup  that  was  used  is  shown  in  figure  5.3.     It  consists  of  two 
wideband  antenna  horns  set  as  a  transmitting  receiving  pair,  an  impulse  generator  feeding 
the  transmitting  horn,  and  a  wideband  sampling  oscilloscope  connected  at  the  output  port  of 
the  receiving  antenna  to  display  the  received  waveforms.     A  layer  of  absorbing  material  was 
placed  such  that  it  would  prevent  reflections  from  the  floor.     When  reflections  were 
desired,  the  absorbing  material  was  replaced  by  a  large  metallic  sheet. 

*Figure  5.2  may  be  compared  with  figure  3.10  to  show  cepstrum  spread  as  compared  to 
clusters . 
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The  antenna  horns  [26]  are  identical  in  their  physical  structure  which  is  shown  in 
figure  5.4.     The  horn  is  designed  in  a  symmetrical  fashion  around  a  center  ground  plane. 
The  antenna  presents  a  50  ohm  driving  point  impedance  which  minimizes  reflections  in  the  50 
ohm  transmission  lines  associated  with  the  transmitter  (impulse  generator)  and  receiver 
(oscilloscope).  , 

The  impulse  generator  has  a  50  ohm  source  impedance  and  generates  the  pulse  shown  in 
figure  5.5.  The  oscilloscope  sampling  head  terminates  the  receiving  transmission  line  in 
50  ohm. 

~  '         '  "      -  '  ■•    ■  5.4    The  Application 

The  output  signal  of  the  receiving  horn  of  the  setup  shown  in  figure  5.3,  with  the 
absorbing  material  in  place,  was  acquired  and  is  shown  in  figure  5.6.     The  figure  shows 
some  reflection  and  scattering  components  not  overlapping  the  main  signal.  Therefore, 
ordinary  time  windowing  can  be  used  to  remove  these  unwanted  components.     However,  the 
homomorphic  deconvolution  transform  will  be  applied  for  the  purpose  of  demonstrating  the 
effect  of  the  separating  filter  on  the  separated  output  signal. 

Applying  the  homomorphic  deconvolution  transform  on  this  signal,  the  cepstrum  of 
figure  5.7  is  obtained,  which,  as  was  discussed  earlier,  contains  a  spread  of  contributions 
rather  than  discrete  clusters.     The  separation  process  was  tried  twice  on  this  cepstrum 
•using  different  filters  having  the  half-widths  A  and  B,  respectively,  as  indicated  in 
figure  5.7  (the  filters  were  symmetric  about  the  cepstrum  origin).     The  result  of  the  first 
separation  at  position  A  is  given  in  figure  5.8,  which  may  be  compared  with  figure  5.9 
where  the  reflections  were  time-windowed  out  at  the  indicated  point  on  the  figure.     The  two 
results  very  closely  match,  while  the  separation  at  B,  shown  in  figure  5.10,  matches  the 
time-windowed  signal  of  figure  5.11.     From  these  figures  it  can  be  concluded  that  the 
narrower  the  window  in  the  cepstrum  the  narrower  the  window  in  the  corresponding  time 
domain.     Also,  a  too  narrow  window  in  the  cepstrum  will  start  distorting  the  main  signal  as 
is  clear  from  separation  B  where  some  distortion  showed  up  just  preceding  the  main  signal, 
figure  5.10b. 

After  experimenting  with  the  above  separation  trials,  it  is  interesting  now  to  try  a 
case  where  the  reflections  overlap  the  main  signal.     This  was  achieved  by  introducing  a 
large  metallic  plate  in  place  of  the  absorbing  material  in  the  setup  of  figure  5.3.  The 
observed  waveform  is  given  in  figure  5.12  where  it  is  possible  to  see  the  reflections  from 
the  metallic  plate  overlapping  the  main  signal.     The  corresponding  cepstrum,  obtained  by 
the  application  of  the  homomorphic  transform,  is  shown  in  figure  5.13.     Two  separations 
were  used  again  at  the  windows  indicated  in  the  figure.     The  results  of  the  two  separations 
are  given  in  figures  5.14  and  5.15.    The  two  separation  results  are  slightly  different  from 
each  other  and  from  the  clean  signal  of  either  figure  5.9  or  5.11,  but  it  is  important  to 
notice  that  the  introduced  overlapping  reflections  (by  the  metallic  plate)  were  successfully 
removed  through  the  use  of  the  homomorphic  deconvolution  transform. 
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5.5    Summary  and  Conclusions 

The  present  chapter  began  with  an  introduction  about  the  role  that  the  time  domain 
measurement  techniques  can  play  in  solving  electromagnetic  radiation  problems.     Then  the 
distributed  nature  of  the  expected  cepstrum  in  such  problems  was  discussed  and  compared 
with  previously  studied  discrete  cases .     Two  examples  of  applying  the  homomorphic  deconvolu- 
tion  transform  to  an  antenna  system  were  presented,  and  the  transform  proved  to  be  very 
successful  in  the  removal  of  the  unwanted  reflections  and  scattering  components  either 
overlapping  or  not  overlapping  the  main  signal  containing  the  desired  signal  information. 
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DIELECTRIC  FILLING 


Figure  5.4.   The  Wideband  Antenna  Horn. 

(all  dimension  in  cm) . 
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6.     TRANSFER  NETWORK  IDENTIFICATION  (MODELING)  OF  FEED-THROUGH  SAMPLING-HEAD 


-■  6.1  Introduction 

In  chapters  3  and  4  it  was  noted  that  a  given  oscilloscope  recording  a  time  domain 
signal  always  presents  a  recorded  waveform  different  from  the  signal  waveform,  which  is  due 
to  the  limited  oscilloscope  bandwidth.     The  recorded  signal  is  the  convolution  of  the  signal 
waveform  with  that  of  the  oscilloscope  impulse  response,  figure  6.1.     Also,  it  was  shown  in 
these  chapters  that  the  oscilloscope  bandwidth  limitation  did  not  prevent  the  application  of 
the  homomorphic  transformation  to  solve  problems  in  which  the  oscilloscope  transfer  function 
was  eliminated  by  insertion  (ratio)  measurements.     However,  if  it  is  necessary  to  determine 
the  true  waveform  at  the  input  of  the  oscilloscope,   then  a  model  for  the  transfer  properties 
of  the  oscilloscope  is  required. 

In  time  domain  spectroscopy,  the  overall  bandwidth  is  determined  by  the  product  of  the 
bandwidths  of  the  signal  source  and  the  oscilloscope.     In  order  to  analyze  the  relative 
bandwidth  contributions,  it  is  necessary  to  determine  a  model  for  the  oscilloscope. 

For  the  picosecond  time  domain  signals  which  are  being  dealt  with  in  this  work,  the 
wideband  sampling  oscilloscope  is  used  for  observing  and  recording  such  signals.     The  one 
kind  of  sampling  head  that  is  most  commonly  used  in  time  domain  ref lectometry  measurements 
is  the  feed-through  sampling-head  because  of  the  convenience  of  usage  in  observing  incident 
and  reflected  time  domain  signals  on  the  signal  line,  figure  6.2. 

The  aim  of  this  chapter  is  to  develop  a  network  model  based  upon  the  physical  properties 
of  the  sampling-head.     Nominally,  the  manufacturer  claims  a  28  picosecond  (10%  to  90%) 
transition  time.     The  network  model  is  to  be  used  to  predict  both  the  oscilloscope's  time 
and  frequency  domain  responses.     The  predicted  model  responses  are  to  be  used  to  deconvolve 
the  observed  signals  in  order  to  deduce  the  true  waveform  at  the  oscilloscope  input  port. 
Previous  workers  have  attempted  to  use  heuristically  based  mathematical  models  which 
considered  only  an  idealized  sampling  process  [27]. 

6.2    Theory  and  Construction  of  Sampler 

Previous  publications   [28,29]  went  in  depth  into  the  theory  and  construction  of  the 
feed-through  samplers;  only  a  brief  coverage  will  be  presented  here. 

Figure  6.3  shows  the  balanced  two-diode  sampling  circuit;   this  circuit  is  the  basic 
structure  in  the  feed-through  sampler.     The  circuit  works  in  the  following  manner.  When 
the  signal  on  the  feed-through  line  is  zero,  application  of  a  balanced  strobing  (sampling 
command)  pulse  on  the  two  diodes  will  result  in  equal  conduction  in  the  diodes;  consequently, 
an  equal  charge  storage  results  on  the  sampling  capacitors  C^  which  in  turn  yields  a  zero 
output.     But  if  the  signal  on  the  line  has  a  nonzero  value,  nonequal  diode  conduction  will 
result  in  a  nonequal  charge  storage  on  the  sampling  capacitors  and  an  output  different  from 
zero.     The  output  polarity  and  magnitude  are  related  to  the  input  polarity  and  magnitude. 
The  function  of  the  feedback  is  to  readjust  the  diodes'  reverse  bias  after  each  sampling 
such  that  the  next  sampling  process  only  detects  the  change  (the  difference)  in  the  input 
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signal.     The  feedback  serves  two  purposes:     (1)  it  increases  the  dynamic  range  for  the 
allowable  input  voltage  to  the  sampler,  and  (2)  it  maintains  the  operating  region  of  the 
diode  almost  constant  during  each  successive  sample. 

Since  the  output  of  the  sampling  circuit  is  a  measure  of  the  differences  or  level 
changes  in  the  input  signal,  an  integrating  circuit  is  needed  to  give  the  required  measure 
of  the  input  signal.     Figure  6.4  shows  the  dynamics  of  the  sampling  network  assuming  a  ramp 
input  signal  and  a  triangular  sampling  pulse. 

The  detailed  view  of  the  physical  realization  of  the  two-diode  feed-through  sampler  is 
shown  in  figure  6.5.     In  the  figure  it  is  seen  that  the  feed-through  transmission  line  is 
tapered  or  stepped-down  on  both  the  input  and  output  sides;  this  is  successively  done  through 
a  three  section  compensated  stepped-taper  maintaining  a  50  ohm  impedance  along  the  line 
(details  shown  in  figure  6.8).     The  sampling  diodes  are  put  in  a  perpendicular  fashion  to 
the  transmission  line  touching  the  inner  conductor  from  both  sides .     In  this  region  the 
inner  conductor  is  a  thin  curved  piece  of  metal  that  joins  the  two  tapered  sections  together. 
The  outer  conductor  forms  a  biconical  cavity,  which  acts  as  a  shorted  biconical  transmission 
line  which  shapes  the  strobe  pulse,  as  shown  in  figure  6.6.     The  biconical  cavity  is  dielec- 
tric filled  to  mask  the  packaging  capacitance  of  the  sampling  diodes  to  lessen  the  effect  of 
the  diode  on  the  50  ohm  characteristic  impedance  throughout  the  cavity. 

The  sampling  capacitors,  2  pf  each,  are  built  in  the  sampling  diodes;  figure  6.7  shows 
the  details  of  a  diode  pack. 

6.3    Modeling  of  Physical  Structure 

Precise  measurements  were  performed  on  the  physical  dimensions  of  a  disassembled 
sampling-head;  the  results  are  given  in  figure  6.8.     Each  of  the  step  line  discontinuities 
can  be  modeled  by  an  equivalent  capacitor  [30,31];  the  value  of  this  capacitance  is  a  func- 
tion of  the  inner  and  outer  diameters  of  the  discontinuity.     The  discontinuity  dimensions 
and  capacitances  are  tabulated  in  Table  6.1.     Figure  6.9  shows  the  equivalent  circuit  for 
the  tapered  structures  complete  with  the  characteristic  impedances  of  line  portions  in 
between. 

6.4    Biconical  Cavity 

The  details  of  the  cavity  have  been  shown  in  figures  6.5,  6.6,  and  6.8.     The  cavity 
parameters  of  interest  are  (1)  the  characteristic  impedance  as  seen  by  the  signal  along  the 
feed-through  line,   (2)  the  characteristic  impedance  of  the  biconical  transmission  line  as 
seen  by  the  strobing  pulse,    (3)  the  travel  time  out  and  back  along  the  shorted  biconical 
line  which  determines  the  shape  and  duration  of  the  sampling  pulse  applied  to  the  two-diode 
sampling  gate,  figure  6.6,  and  (4)  the  characteristic  impedance  of  the  wire  line  carrying 
the  strobing  pulse  down  to  the  cavity. 
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Table  6.1. 


Discontinuity  Dimensions    (mm)  Capacitances  (pf) 


a  " 

b 

c 

d 

C  . 

m 

C 

out 

6  . 

9977 

3 

.0429 

4. 

9835 

2 

.  1636 

.00607 

.00884 

4. 

9835 

2 

.  1636 

2  . 

8575 

1 

.  2413 

.00803 

.01016 

2. 

8575 

1 

.2412 

2. 

0320 

0 

.8890 

.00246 

.00349 

2. 

0320 

0 

.8890 

0 

.6604 

.00170 

2. 

0447 

0 

.8814 

0 

.6604 

.00170 

2. 

8448 

1 

.  2446 

2. 

0447 

0 

.8814 

.00246 

.00349 

4. 

9962 

2 

.1654 

2. 

8448 

1 

.  2446 

.00803 

.01016 

7  . 

0104 

3 

.0218 

4. 

9962 

2 

.  1654 

.00607 

.00884 

7  . 

0104 

3 

.0218 

2 

.0244 

.00780 

7  . 

0104 

3 

.0231 

2 

.0244 

.00780 
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To  find  (1),   (2),  and  (3)  above,  two  TDR  measurements  were  performed  on  the  cavity  with 
and  without  the  dielectric  filling,  with  the  diodes  removed  in  both  cases.     The  results  are 
shown  in  figure  6.10.     Figure  6.10a  clearly  shows  that  the  feed-through  line  is  within 
±  1  ohm  from  the  50  ohm  characteristic  impedance;  i.e.,  the  reflected  voltage  is  no  greater 
than  ±  5  millivolts.     Also,  it  is  seen  that  there  is  an  inductive  effect  in  the  central 
portion  of  the  cavity  due  to  the  metallic  curved  strip  that  connects  the  tips  of  the  two 
inner  conductor  tapered  sections.     From  figure  6.10b,   it  is  found  that  the  absence  of  the 
dielectric  filling  caused  a  change  in  the  impedance  level  that  corresponds  to  a  reflection 
coefficient  of  p  =  0.196,  but  since  p  =  (v^-1) / (/k+1) ,  where  k  is  the  dielectric  constant,  k 
is  obtained  from        =  (l+p)/(l-p)  =  1.49  which  gives 

k  =  2.21.  (6.1) 

The  value  of  the  dielectric  constant  k  is  to  be  used  to  calculate  the  travel  time  , 

be 

figure  6.6.  From 

=  ^  d^  /c  (6.2) 
be  be 

d^^  is  the  corresponding  travel  distance  which  is  equal  to  4.48  mm  as  given  in  figure  6.8, 
and  c  is  the  speed  of  light;  therefore 

1 1 

=  1.49  X  4.48/(3x10     )  =  22.2  ps.  (6.3) 

DC 

Tne  third  cavity  parameter  to  be  determined  is  the  characteristic  impedance  Z^^  of  the 
"shorted"  biconical  line  which  can  be  calculated  using  the  relation  [32], 

Z^^  =  (l20  in  cot  (a/2)  J /vie  ohms  (6.4) 

where  i  is  the  bicone  half  angle  which  is  equal  to  38° 22'   as  given  in  figure  6.8.  The 
computed  value  of  Z^^  is 

Z^^  =  85.2  ohms.  (6.5) 

Turning  now  to  the  wire  line  which  carries  the  strobing  pulse  to  the  biconical  cavity, 

the  characteristic  impedance  Z    of  the  wire  line  is  to  be  estimated.     The  inner  conductor  of 
^  w 

the  line  has  an  outer  diameter  r^  =  0.127  mm,  coated  with  varnish  to  a  diameter        =  0.33  mm, 

while  the  inner  diameter  of  the  outside  conductor  is  r^  =  0.432  mm.  An  estimate  for  the 

value  of  the  dielectric  constant  k    of  the  used  varnish  is  k    =2.     The  impedance  Z  can 

V  V  w 

be  written  as  the  sum  of  the  two  impedances  Z  ,   and  Z  ^. 

^  wl  w2 

Z    =  Z  ,  +  Z  -  =  — -  In  —  +60  in  — 

W  wl  W2  /r- 

/k         r^  r„ 
V  1  2 

60  „     0.33    ^        .  0.43 
=  —  In    +  60  j^n   

/2        0.127  0.33 
=  57  ohms.  (6.6) 


6.5    Sampling  Diodes 

The  sampling  diodes  that  are  being  used  In  the  feed-through  sampler  are  hot  carrier 
(Schottky)  diodes;  the  electrical  model  of  such  diodes   [33,34,35]  Is  shovm  In  figure  6.11, 
in  which   ,..  ,  ;  ,  -  . 


I 

s 

R. 

J 
C. 

2 

R 

s 

^d 
C 

P 


the  diode  (reverse)  saturation  current, 

the  junction's  dynamic  resistance, 

the  junction  capacitance, 

the  spreading  and  contact  resistance, 

the  lead  inductance,  and 

the  package  capacitance, 


wnere 


R.  =  1/ 


31. 


'•3V, 


(6.7) 


and 


J 


exp(v./(f^)l   -  1 


:fl-.. 


(6.8) 


In  the  following,  the  measurements  required  for  the  evaluation  of  the  above  diode 
parameters  are  divided  into  two  types:     (1)  static  to  evaluate        and        and  (2)  dynamic  to 

evaluate  C   ,  L  , ,  and  C.. 

P      d  J 

6.5.1    Static  measurements  (I    and  R  ) 

s  s 

The  circuit  of  figure  6.12  was  used  to  measure  the  diode's  dc  characteristics  by  taking 

values  of  the  dc  current  I,  for  various  values  of  dc  diode  voltage  V,;  the  measured  values 

Q  d 

are  plotted  in  figure  6.13.     In  this  dc  measurement,  the  diode  currents  and  voltages  are 


^j  =  ^d' 


(6.9) 


V_,  -  I.R  . 
d        J  s 


(6.10) 


In  the  forward  bias  case  where  I.  >>  I  ,  the  diode  I-V  characteristic,  eq.   (6.8),  can  be 

J  ^ 

approximated  as 


I.  =  I  exp(V./<))) 


(6.11) 


£n  I.  =  Jin  I    +  V./(j);       I.  >>  I 
J  s        J  ^  J  s 


(6.12) 


where 

(p  =  kT/e. 

Substituting  eq.   (6.9)  in  eq.   (6.12)  gives 


Jin  I,  =  Jin  I    +  (V,-I,R  )/(});      I,  >>  I  . 
d  s  dds  d  s 


(6.13) 
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The  experimental  data  of  figure  6.13  were  fitted  to  eq.   (6.13)  above  yielding  the  following 


values  for  R    and  I  : 
s  s 

=  19  ohms  (6.14) 
=  2.76  X  10        amps.  (6.15) 

6.5.2    Dynamic  measurements  (C  ,  L,,  and  C.) 

^  p      d  J 

To  determine  the  sampling  diodes'  dynamic  parameters  (C^,  L^,  and  C_. ) ,  the  TDR  tech- 
nique was  used;  the  setup  for  the  experiment  is  shown  in  figure  6.14.     The  test  signal  was 
the  tunnel  diode  pulse  given  in  figure  6.15.     Because  of  the  limited  resolution  due  to  the 
nonzero  risetime  of  the  test  signal,  the  performed  TDR  measurements  were  not  able  to  directly 
resolve  the  contributions  of  the  sampling  diodes  from  those  of  the  central  conductor  inductanc 
L^;  i.e.,  the  TDR  views  the  sampling  diodes  and  the  central  inductance  as  if  they  were  at 
the  same  physical  location.     This  lumping  can  be  tolerated  (and  later  resolved)  since  the 
metallic  strip  extends  for  about  1  mm  on  both  sides  of  the  diodes  in  a  symmetrical  fashion. 

This  symmetry  allows  L    to  be  divided  equally  into  two  equal  parts,        and  L„  on  the  left 

C  L  K 

and  right  of  the  diodes,  respectively,  figure  6.16. 

To  overcome  the  TDR  lumping  of        and  the  diodes'  contributions,  two  TDR  measurements 
were  performed;  one  with  both  diodes  in  place,  and  the  other  one  with  only  one  of  the  diodes 
in  place,  figures  6.17a  and  6.17b.     Denoting  the  maximum  voltage  drop-peaks  as  e^  and 
respectively,  the  TDR  voltages  are 

^  =  2e^  +  (6.16) 

62  =        +  (6.17) 

where  e^  and  e^  are  the  contributions  in  the  TDR  signal  corresponding  to        and  diode, 
respectively.    Solving  eqs.  (6.16)  and  (6.17)  for  e    and  e,  yields 

Li  Q 

e^^  =  2e2  -  e^  (6.18) 
=        -  e^.  (6.19) 


Substituting  for  e^  and  e^  from  figure  6.17  gives 


e^  =  -9.9  mV  (drop-peak)  =  9.9  mV  peak  (6.20) 

e^  =  19.7  mV.  (6.21) 

The  value  of  eq.   (6.20)  can  be  used  to  evaluate  the  inductance        using  the  relation  [36] 
2Z 

L  =  —  e-  (6.22) 

m      L  ^  ' 

where  m  is  the  maximum  slope  of  the  test  pulse,  figure  6.15.  Therefore, 

=  ^^^-^    ^  9.9  X  10"^  =  150  ph  (6.23) 
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and 

L^=  75  ph.  (6.24) 

To  measure  the  package  capacitance        of  the  diode,  a  third  TDR  measurement,  figure 
6.17c,  was  done  with  one  working  diode  and  an  open  diode;  the  open  diode  contributes  only 
Denoting  the  voltage  drop-peak  in  this  measurement  e^,  the  TDR  voltage  relation  is 

^3  =  ^d  +  %  +  .  (6.25) 

where  e    is  the  contribution  of  the  open  diode,  or  in  other  words,  the  C  contribution. 

P  P 
Consequently,  ., 

%  =  ^3  ■  ^d  "  '  (6.26) 

Upon  substituting  for  e^  from  figure  6.17c,  and  for  e^^  and  e^  using  eqs.   (6.20)  and  (6.21), 

the  value  of  e     is  given  by 

P 

e    =  20.8  -  19.7  +  9.9  =  11  mV.  (6.27) 

can  then  be  calculated  using  [36]  ■•■n 


2_ 
"o 


which  yields 

2  •--      -  , 

C    =   X  11  =  .067  pf  .  (6.29) 

P      6.58  X  50 

For  the  order  of  magnitudes  involved  in  the  present  TDR  measurement,  computations 

showed  that  the  effect  of  diode  inductance        on  the  voltage  drop-peaks  is  minor;  i.e., 

such  measurements  are  not  able  to  yield  the  value  L^.     However,  previous  measurements 

[37]  that  were  done  on  the  same  particular  sampler  gave  a  complex  value  14  +  j5  ohms  at 

10  Griz  for  the  impedance  Z^/2  of  the  sampling  gate  with  a  diode  biasing  forward  current  of 

5  ma.     This  information  can  be  used  to  calculate  the  value  of  L,  as  follows:     the  measured 

d 

impedance  at  the  sampling  gate  is  in  fact  the  impedance  of  the  parallel  combination  of  the 
two  sampling  diodes,  each  one  in  series  with  its  sampling  capacitor  C^;  this  combination  is 
shown  in  figure  6.18.     At  5  ma  diode  forward  current,  the  junction  resistance  is 


R.  =  1/ 

J 


81. 


3V. 
J 


<^/l^  =  5.2  n  (6.30) 


R./2  =  2.6  n.  (6.31) 
J 

Since  C.  is  expected  to  be  lower  than  0.1  pf,  the  reactance  of  2C^  at  10  GHz  is  higher  than 
80        which  is  much  greater  than  the  reactance  of  2C^  can  be  neglected  in  this 
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particular  case.    Hence,  the  impedances  can  be  expressed  as 
Z  (l/ja)2C  )   •   (R  /2+R./2+ja)L,/2) 

_s.  1_  +   2  §  a  3         .  (6,32) 

2        jw2C        (l/j(jj2C  )  +  (R  /2+R./2+ja)L,/2) 

S  P  SJ  Q 

Since  C  ,  R  ,  and  R.  are  known,  eq.   (6.32)  above  can  be  solved  to  give 
p      s'  J 

=  300  ph.  (6.33) 

Now,  the  only  component  that  is  left  to  be  determined  is       ;  and  since  it  is  known  to 
be  a  function  of  the  junction  voltage,  a  set  of  TDR  measurements,  with  both  diodes  in  place, 
was    performed  at  different  biasing  voltages.     The  biasing  circuit  was  that  shown  in  figure 
6.19,  and  in  figure  6.20  only  two  of  the  various  TDR  waveforms  are  shown;  however,  all  the 
results  are  listed  in  table  6.2.     Denoting  the  voltage  drop-peak  in  these  TDR  runs  by  e(v), 
where  v  is  the  diode  bias,  the  TDR  voltage  is 

e(v)  =  2e    +  2e.(v)  +  e^  (6.34) 
P         J  ^ 

from  which  e.(v),  the  contribution  due  to  C . ,  can  be  obtained  as 
J  J 

e^(v)  =  [e(v)  -        -  2e^]/2.  (6.35) 

The  corresponding      (v)  is 

C.(v)  =  -f-  e  (v).  (6.36) 

These  steps  are  shown  also  in  table  6.2,  and  the  calculated  C.(v)  is  plotted  in  figure  6.21. 
By  evaluating      (v) ,  this  phase  of  modeling  the  sampling  diodes  is  concluded. 
In  summary  the  obtained  values  for  the  diode  model,  figure  6.11,  are 

C    =  0.067  pf 
P 

=  300  ph 

R    =  19  n 
s 

Ig  =  2.76  X  lo"-*--"-  amp 

(v)  is  given  in  figure  6.21.  (6.37) 
Also,  the  value  of  the  central  conductor    inductances  has  been  determined  as 

L-  =  L„  =  75  ph.  (6.38) 

6.6    The  Sampling  Pulse 

The  experimental  setup  of  figure  6.22  was  used  to  observe  the  strobe  pulse  waveform; 
only  the  leading  edge  of  the  pulse  will  be  considered  since  it  is  the  only  part  of  the 
waveform  that  forms  the  sampling  pulse,  figure  6.6.     The  observed  signal  is  shown  in 
figure  6.23  which  shows  an  amplitude  of  9.58  volts  and  a  10%  to  90%  risetime        of  53.5  ps; 
both  of  these  values  have  to  be  corrected  as  explained  in  the  following  discussion. 
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Table  6.2. 


V  V 

e(v)  mV 

e  ,  (v)  mV 
J 

C.(v)  pf 

3 

0.0 

26.86 

7  .40 

.045 

-0.25 

24  .  90 

6  .42 

.039 

-0.50 

23.58 

5  .76 

.035 

-0.75 

23.26 

5.60 

.034 

-1.00  , 

2  2.60 

5.27 

.032 

-1.25 

21.94 

4.94 

.030 

-1.50 

21.60 

4.77 

.029 

-1.75 

21.28 

4.  61 

.028 

-2.00 

20.62 

4.28 

.026 

-2.25 

20  .63 

4.28 

.026 

-2.50 

20.28 

4.11 

.025 

-2.75 

19.96 

3.  95 

.  024 

-3.00 

19.62 

3.78 

.023 
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The  characteristic  impedances,  Z        Z  „,  and  Z  „,  of  the  stripline  (on  which  the  strobe 

ol      o2  o3 

pulse  was  measured) ,  the  (wire)  microminiature  coaxial  line  connecting  the  stripline  to  the 
biconical  cavity,  and  the  biconical  line  are  50,  57,  and  85.2  ohms,  respectively,  by  virtue 
of  eqs.   (6.3)  and  (6.4).     These  unmatched  line  impedances  cause  a  voltage  amplification 
factor  A  given  by 

A  =  (l+p^)(l+P2)  (6.39) 
where  the  reflection  coefficients        and        are  determined  as 

^o2      ^ol      57  -  50  ,^  , 

=  z  ,  +  z  ,  =  57T-50  =  -o^^  ^^-"^^^ 

02  ol 
Z      -  Z 

03  o2      85.2  -  57      „  mo  l^\ 

'2  =  Z  ,  +  Z  ,  =  85.2  +  57  =  ^'^^^  ^^"^'^ 
o3  o2 

from  which  A  is  calculated  as 

A  =  (1+.065)   •   (1+0.198)  =  1.276.  (6.42) 

The  amplitude  of  the  observed  strobe  pulse  vs.  time  has  to  be  multiplied  by  the  factor  1.2  76 
to  yield  the  corresponding  amplitude  at  the  input  of  the  biconical  shorted  line. 

The  risetime  correction  has  to  take  into  consideration  the  attenuators,  adaptors,  and 
measuring  sampling-head  responses.     Experiments   [38]  have  shown  that  elements  have  the 
following  nominal  risetimes: 

20  ps  for  10  dB  attenuator, 
18  ps  for  20  dB  attenuator, 
10  ps  for  each  adaptor,  and 
20  ps  for  measuring  sampling-head. 

Using  the  risetime  square  rule  for  the  cascaded  setup  of  figure  6.22,  the  equivalent  rise- 
time  T    of  this  combination  is  nominally 
e 


T    =  /202  +  2  X  io2  +  202  +  18^  =  36.4  ps.  (6.43) 
e 

Using  the  square  rule  again,  the  estimated  risetime        of  the  strobe  pulse  is 


X    =  /t2    _  x2    =  /53.52  -  36. 42  =  39.2  ps,  (6.44) 
s  m  e 

An  assumed  response  having  36.4  ps  risetime  could  be  used  to  deconvolve  the  observed  strobe 
pulse  to  yield  an  estimate  of  the  leading  edge  of  the  actual  strobe  pulse  waveform.  Another 
method,  which  was  found  to  be  as  good,  is  to  scale  the  time  variable  in  the  observed  wave- 
form by  the  ratio  t  /i     to  yield  an  estimate  for  the  actual  waveform.     The  two  methods  gave 
s  m 

values  which  were  within  5%  of  each  other. 

The  corrected  leading  edge  waveform  for  the  strobe  pulse  at  the  input  port  of  the 
shorted  biconical  line  is  given  in  figure  6.24. 
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:  6.7     Structure  Step  Response 

In  sections  6.2  through  6.6  all  the  sampling-head  components  were  modeled.  The  com- 
plete model  (equivalent  circuit)  would  be  the  same  as  that  shown  in  figure  6.9  except  for 
the  cavity  box  to  be  replaced  by  its  model,  figure  6.25. 

To  study  the  step  response  of  the  sampler,  the  computer  program  TDNA  given  in  appendix  C 
is  to  be  used;  a  voltage  step  is  to  be  applied  to  the  input  port,  and  the  sampled  output  is 
to  be  integrated  to  yield  the  step  response.     As  seen  from  figures  6.9  and  6.25,  the  number 
of  elements  in  the  model  is  quite  large;  consequently,   the  amount  of  computer  memory  re- 
quired to  consider  this  model  was  beyond  the  capacity  of  a  large  computer.     The  network 
(model)  of  figures  6.9  and  6.25  was  divided  into  some  independent  networks,  in  which  case 
the  step  response  can  be  obtained  in  two  steps:     the  structure  step  response  and  the  sampling 
process.     Figure  6.26  shows  the  simplified  network  considered  in  the  first  step.     Here  the 
short  lines  of  figure  6.9  and  the  discontinuity  capacitors  on  each  side  are  replaced  by  one 
lumped  capacitor,  and  two  lines  on  both  sides  were  matched  in  impedance  to  their  neighboring 
lines  in  impedance.     In  this  way  the  number  of  lines  and  capacitors  are  reduced  to  about 
one  half.     Also,  in  the  figure  the  sampling  signal  and  feedback  voltage  are  ignored;  the 
diode  bias  is  set  to  a  constant  negative  value        =  -1.63  volts,  at  which  the  actual  sampling 
diodes  normally  operate,  as  was  experimentally  measured. 

The  network  parameters  shown  in  figure  6.26  were  given  to  the  computer  as  input  data 
for  the  program  TDNA.     The  obtained  time  domain  signal  Vg  at  the  node  8  is  shown  in  figure 
6.27. 


6.8    Sampling  Process 

The  second,  and  final,  step  in  obtaining  the  sampling-head  step  response  is  to  sample 

the  resultant  voltage  waveform  of  node  8  using  the  computer  program  SMPL  given  in  appendix  C. 

The  voltage  signal  V„  is  to  be  applied  to  node  17  in  figure  6.28,  each  time  the  balanced 
o 

strobing  signal  is  applied  to  nodes  2  and  16.     The  reverse  bias  of  the  diodes        is  kept 

at  -1.63  volts,  the  normal  operation  value.     The  difference  between  the  resultant  voltage 

increments  on  the  sampling  capacitors  C  ,  and  C  „  is  a  measure  of  the  output  voltage  V 

si  s2  ^  ^  out 

The  sampling  efficiency  n   [27,28]  is  defined  as 
m 

As  mentioned  earlier  in  section  6.2,  the  feedback  voltage  readjusts  the  diodes'  reverse 
bias  after  each  sample  such  that  the  next  sampling  process  only  detects  the  change  in  the 
input  signal;  consequently,  at  the  beginning  of  the  sampling  interval,  the  net  diode  bias 
is  a  constant  -1.63  volts.     Accordingly,  and  referring  to  figure  6.28,  the  value  of  V 


fb 


for  a  new  sampling  process  should  be  equal  to  V^^  of  the  previous  sampling; 


V^^  =  V.   .  (6.46) 
fb  in 
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Figure  5.29  shows  a  block  diagran  of  the  sampling  network  and  its  feedback  from  the 
figure,  and  using  eqs .   (6.45)  and  (6.46)  it  can  be  written  that 

fb  out 

=    SG  n  V. 

m 

=  V.   ;  (6.47) 

in 

therefore, 

BGn  =  1  (6.48) 

or 

-,       V.  ' 

BG  =  ^  =  ^.  (6.49) 

out 

To  determine  the  factor  5G,  a  preliminary  run  for  the  computer  program  SMPL  was  done  using 
a  small  voltage  step  at  node  17  as  an  input,  and  upon  calculating  its  ratio  to  the  output 
voltage,   the  obtained  values  were 

V.     =  .01  volt 
m 

V        =  .0026907  volts 
out 

=  .0026907  =  (6-50) 

out 

The  determination  of  the  factor  £G  allows  the  calculation  of  the  feedback  voltage  V„ 

ZD 

after  each  sample  using  eq .  (6.47). 

The  program  SMPL  was  run  using  the  voltage  Vg  as  the  input  signal  to  the  sampling 
network,  and  the  input  signal  was  sampled  at  1  ps  time  intervals.     The  obtained  output  was 
normalized  to  yield  the  required  sampling-head  step  response,  figure  6.30.     The  step  response 
shows  a  30.5  ps  transition  time. 

Figures  6.31  through  5.34  show  the  voltage  waveforms  of  the  sampling  process  for  one 
sampling  interval.     Figure  5.31  shows  the  time  domain  waveforzi  of  the  resultant  sampling 
pulse  ^2-±  (voltage  difference  between  nodes  3  and  1  in  figure  6.28).     The  voltage  across 
one  of  the  diodes        ^  is  shown  in  figure  6.32,  while  the  voltage  across  one  of  the  sampling 
capacitors        ^  is  shown  in  figure  6.33.     The  difference  between  the  two  sampling  capacitor 
voltages  ^2.2-1  ~  ^3-6  ''"^  sampled-data  output,  figure  5.3^.     It  is  seen  that  the  resultant 

shape  of  the  sampling  pulse  is  somewhat  different  from  being  triangular,  and  the  sampling 
time        (the  time  of  diode  conduction)  is  about  21  ps . 

In  summary,  it  has  been  shown  that  the  waveform  of  the  sampling  process  or  "sampling 
pulse"  is  derived  from  the  strobe  pulse  step  like  transition  through  shaping  by  the  biconical 
cavity  into  a  quasi-triangular  pulse  which  is  in  turn  shaped  by  the  nonlinear  and  band- 
limiting  properties  of  the  sampling  diode  pair. 
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■  6.9    The  Effect  of  the  Diode  Bias  on  the  Sampling-Head  Step  Response 

For  ideal  sampling,  the  sampling  pulses  should  be  infinitesimal  in  duration;  a  dirac- 
delta  function  is  the  ideal  shape.     The  actual  sampling  pulse,  figure  6.31,  does  have 
nonzero  duration;  consequently,   the  step  response  of  the  sampling-head  is  expected  to  have 
a  nonideal  (nonzero)  transition  time.     From  the  theory  of  sampling  processes  it  is  clear 
that  the  greater  the  duration  of  the  sampling  pulse,  the  slower  will  be  the  step  response 
transition  time. 

There  is  one  way  for  the  operator  to  control  the  sampling  pulse  duration  and,  consequently, 
the  step  response;  this  is  accomplished  by  controlling  the  diodes  reverse  bias,  thus  varying 
the  level  at  which  the  diodes  conduct  [39],  figure  6.35.     The  larger  the  diode  reverse 
bias,  the  smaller  will  be  the  sampling  time;  consequently,  the  faster  will  be  the  response 
of  the  sampling-head,  and  vice  versa.     However,  there  are  limits  to  bias  variations;  very 
low  or  very  high  reverse  bias  will  cause  imperfect  sampling.     Figure  6.36  shows  the  observed 
waveforms  of  a  tunnel  diode  step  using  the  feed-through  sampling-head  for  three  values  of 
diode  reverse  bias  corresponding  to  (a)  bias  for  optimum  sampling,   (b)  low  reverse  bias, 
and  (c)  high  reverse  bias. 

Additional  tunnel  diode  waveforms  were  measured  for  various  values  of  the  diode  reverse 

bias.     The  10%-90%  transition  time  t    was  calculated  for  each  waveform.     The  obtained  t 

r  r 

was  plotted  vs.  reverse  bias  as  in  figure  6.37.     Since  the  tunnel  diode  transition  time  is 
constant,  the  observed  variations  are  due  to  the  sampling-head  response  variations.  The 
values  of  the  sampling  time        for  various  values  of  reverse  bias  were  computed  (using  the 
program  SMPL)  and  also  plotted  in  figure  6.37.     The  figure  agrees  with  the  fact  that  the 
smaller  the  sampling  time  the  faster  the  response  of  the  sampling-head,  within  the  limits 
of  effective  feedback  control  on  the  sampling  level  voltage.     For  reverse  bias  values 
higher  than  1.76  volts,  the  response  is  slower  unlike  the  theoretical  expectations.  From 
the  figure  it  is  clear  that  the  fastest  response  for  the  sampling-head  can  be  obtained  by 
adjusting  the  diode  reverse  bias  to  the  value  of  1.76  volts;  the  observed  waveform  at  this 
bias  was  slightly  different  from  the  normal  sampling  case  (1.63  volts  bias  which  results 
from  the  manufacturer's  alignment  procedure  for  proper  operation). 

The  sampling-head  step  response  for  the  fastest  response  bias  (1.76  volts)  was  com- 
puted and  is  given  in  figure  6.38,  where  it  is  seen  that  the  step  response  shows  a  27.5  ps 
transition  time.     Figure  6.39  gives  the  corresponding  impulse  response  (time  derivative). 

6.10    The  Deconvolution  Routine 

The  final  step  in  the  present  chapter  is  to  develop  a  deconvolution  routine  that  uses 
the  sampling-head  impulse  response  given  in  figure  6.39  for  deconvolving  an  observed  i 
signal  in  order  to  deduce  the  true  waveform  at  the  sampling-head  input  port.     In  other 
words,  and  referring  to  figure  6.1,  the  problem  is  to  find  x(t)  knowing  both  y(t)  and  h(t). 

The  simplest  deconvolution  method  known  is  to  calculate  the  ratio  of  the  frequency 
domain  transforms  ¥(6"^^) /H(e~''^)   to  obtain  X(e-''^),   the  frequency  domain  form  of  x(t)  .  This 
method  inherently  produces  a  very  noisy  output,  because  in  the  division  process  the  zeros 
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and  the  poles  of  H(e'^'^)  must  be  exactly  cancelled  out.    Any  errors  in  the  data  representing 
Y(e"''^)  and  H(e'''^)  degrade  the  cancellation  and  introduce  rapid  variations  (noise)  in  the 
resultant  X(e"''^) .     In  order  to  suppress  the  noise,  low  pass  filtering  and/or  pole  removal 
processes  have  to  be  applied,  the  latter  acting  on  poles  introduced  by  incomplete  pole/zero 
cancellations.     The  pole  removal  processes  are  very  sensitive  which  makes  the  use  of  the 
simple  ratio  method  for  deconvolution  rather  unlikely. 

Although  more  complicated  deconvolution  methods  are  available  for  use,  they  all  involve 
successive  approximations  to  an  estimated  input  x(t)  to  minimize  the  error  between  the 
estimated  output  and  the  actual  output  y(t) . 

In  this  section  a  slightly  modified  homomorphic  deconvolution  method  will  be  used 
successfully  to  deconvolve  the  output  y(t)  to  obtain  x(t) .     The  method  is  advantageous  in 
its  simplicity  as  well  as  yielding  a  low  level  of  output  noise  in  x(t) . 

6.10.1    The  homomorphic  deconvolution  method 

In  the  homomorphic  deconvolution  as  was  developed  in  chapter  2,  the  cepstrum  of  a 
convolution  is  the  linear  sum  of  the  cesptras  of  the  individual  components  of  this  con- 
volution.    That  is,  using  the  notation  of  chapter  2  as  applied  to  the  signals  of  figure  6.1, 
it  can  be  written  that 


y(t)  =  x(t)  +  h(t) .  (6.51) 

Knowing  both  y(t)  and  h(t),  the  cepstras  y(t)  and  h(t)  can  be  used  to  determine  x(t) 
through  the  simple  subtraction 

x(t)  =  y(t)  -  h(t).  (6.52) 

The  resultant  x(t)  can  be  processed  (inverse-homomorphically  transformed)  back  to  the  time 
domain  to  yield  the  required  x(t) . 

6.10.2    An  application 

To  show  the  capability  of  the  homomorphic  deconvoluti  on  method  as  presented  above,  the 
method  will  be  applied  to  a  known  convolution.     This  known  convolution  y(t)  was  formed  by 
convolving  the  sampling-head  impulse  response  h(t),  figure  6.39,  and  an  assumed  typical 
step-like  time  domain  signal  x^(t),  figure  6.40.     The  resultant  convolution  y(t)  is  shown 
in  figure  6.41.     The  cesptras    y(t)  and  h(t)  given  in  figures  6.42  and  6.43,  respectively, 
were  obtained  by  applying  the  homomorphic  transformation  to  their  corresponding  time  domain 
forms.     The  difference  x^(t)  =  y(t)  -  h(t)  is  shown  in  figure  6.44.     Applying  the  inverse 
homomorphic  transform  to  x^(t)  yields  x^(t),  the    deconvolved  estimate  for  the  assumed 
input  signal  x  (t).     Comparison  of  figures  6.45  and  6.40  shows  that  x,  (t)  matches  to  a 

3i  Q 

great  extent  the  original  signal  x  (t)  except  for  a  very  small  amount  of  error  oscillations 

3. 

which  are  due  to  sampling  and  computational  errors. 
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6.11    Summary  and  Conclusions 

The  present  chapter  was  devoted  to  the  modeling  of  the  feed-through  sampling-head. 
The  theory  and  construction  of  the  sampling  head  were  reviewed.     The  physical  dimensions  of 
the  structure  discontinuities  were  used  to  formulate  an  equivalent  network  for  the  sampling- 
head  input  and  output  lines.     The  dimensions  of  the  biconical  cavity  together  with  TDR 
measurements  were  used  to  model  the  cavity  part  of  the  sampler.     The  dc  characteristics  of 
the  diode  together  with  TDR  measurements  vs.  sampling  diode  bias  were  used  to  evaluate 
various  diode  model  parameters.     The  strobing  pulse  was  measured,  and  the  data  corrected 
for  both  amplitude  and  risetime. 

All  of  these  sampling-head  parameters  were  simulated  in  the  form  of  a  digital  computer 
program  which  was  run  to  obtain  a  theoretical  estimate  of  the  sampling-head  step  response. 
The  program  was  implemented  in  two  steps  because  of  the  large  number  of  elements  and  infor- 
mation required  to  simulate  the  sampling  system. 

The  step  response  was  found  to  vary  with  the  diode  reverse  bias;  the  higher  the  reverse 
bias  the  faster  was  the  step  response,  within  limits.  The  step  response  at  the  manufacturer 
suggested  bias  and  at  the  bias  for  the  fastest  step  response  was  obtained. 

The  homomorphic  deconvolution  transform  developed  in  chapter  2  was  slightly  modified 
in  order  to  be  used  to  deconvolve  the  sampling-head  response  contribution  from  the  output 
of  the  sampling  oscilloscope.     Deconvolution  of  the  observed  waveform  provides  an  estimate 
of  the  actual  (true)  waveform  at  the  input  port  of  the  oscilloscope. 
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SIGNAL  OF 
INTEREST 


^  OBSERVED 
y(t)  (RECORDED) 
*  SIGNAL 


y(t)  =  x(t)  *  h(t) 
Y(e^^)  =  X(e^^)  •  H(e^^) 

Figure  6.1     Oscilloscope  Response  Effect 
on  Measured  Signals. 
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Figure  6.2.    Feed-Through  Sampling  Oscilloscope 
as  Used  in  TDR  Measurements. 
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Figure  6.3.  Two-Diode  Sampling  Circuit. 

> 
i 


-134- 


DISPLAY 
TIME 


Figure  6.4.    Sampling  Process. 
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Figure  6.6.   Strobe  Pulse  Shaping. 
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Figure  6.7.   Diode  Pack. 
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Figure  6.8.    Sampling  Head  Dimensions 
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Figure  6.9.   Tapered  Line  Equivalent  Circu 
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a.     With  Dielectric. 


With  no  Dielectric. 


gure  6.10.  TDR  on  Cavity  with  no  Diodes. 

Vert.  10  mV/div  (p  =  .04  2  5/div) 
Horiz.   100  ps/div 


-141- 


<!)  =  —  =  .0259  VOLTS  (AT  300° k) 


Figure  6.11.    Diode  Electrical  Model. 
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Figure  6.12.  Circuit  Used  to  Measure  Static 

Characteristics  of  a  Sampling  Diode. 
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Figure  6.13.    Didoe  Static  Characteristics. 
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a.     Horiz.     1  ns/div,  Vert.     5  0  mV/div. 


b.     Horiz.     10  ps/div,  Vert.     50  mV/div. 


Figure  6.15.    TD  Pulse. 
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a.     Two  Diodes 


b.     One  Diode 


Figure  6.17.    To  be  continued  on  next  page. 
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c .     One  Diode  and  an  Open  Diode 


Figure  6.17.    TDR  on  Diodes. 

Vert.   5  mV/div 
Horiz.   100  ps/div 


-148- 


Figure  6.18.  Sampling  Gate  Impedance  Z  /2 . 
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a.     V  =  0.0  Volt. 


b.     V  =     -1.0  Volt. 


Figure  6.20.    To  be  continued  on  next  page. 
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Figure  6.20.   TDR  vs.  Bias. 

Vert.   5  mV/div 
Horiz.   100  ps/div 
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Figure  6.23.   Observed  Strobe  Pulse. 

Vert.     1.581  v/div. 
Horiz.   20  ps/cm. 
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Figure  6.26.  Simplified  Network  for  Structure 
Step  Response. 
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Figure  6.28.    Sampling  Network. 
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Figure  6.29.    Sampling  Efficiency  and  Feedback. 
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a.     Low  Reverse  Bias   (-1.46  Volt). 


b.     Bias  for  Optimum  Sampling   (-1.63  Volt). 


Figure  6.36.    To  be  continued  on  next  page. 
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c.     High  Reverse  Bias   (-1.76  Volt) . 


Figure  6.36.    Effect  of  Diodes'  Reverse 
Bias  on  Sampling. 
Vert.   50  mV/div. 
Horiz.   20  ps/div. 
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7.     SUMMARY  AND  CONCLUSIONS 

This  report  was  devoted  to  the  study  of  the  theory  and  application  of  the  homomorphic 
transformation  to  time  domain  spectroscopy  and  scattering  problems . 

Chapter  1  was  an  introductory  chapter  which  reviewed  the  contents  of  the  work  presented 
in  this  report. 

Chapter  2  started  with  a  review  of  the  convolution  relations  in  both  the  time  and 
frequency  domains  and  the  need  for  deconvolution  (separation)  techniques.     The  general 
class  of  the  homomorphic  systems  was  defined,  and  their  properties  were  stated  showing  how 
they  can  be  used  for  signal  separations.     In  particular,  the  homomorphic  deconvolution 
transform  was  developed.     The  related  problems,  such  as  the  nonanalyticity  of  the  complex 
logarithmic  function  and  the  effect  of  the  linear  phase  component,  were  presented  and 
discussed,  and  solutions  to  such  problems  were  given.     The  computational  realization  of  the 
deconvolution  transformation  was  made  possible  in  the  form  of  digital  computer  subroutines 
HDT  and  HDF. 

Chapter  3  presented  several  classes  of  time  domain  measurement  problems  in  the  field 
of  linear  microwave  network  analysis.     The  homomorphic  deconvolution  transform,  that  was 
developed  in  chapter  2,  was  used  in  deconvoluting  various  convolutions  in  these  problems. 
The  examples  that  were  considered  are : 

(1)  pure  resistive  loading  along  a  transmission  line, 

(2)  commensurate  transmission  line  networks, 

(3)  TDT  type  measurement  and  with  two  reactive  discontinuities,  and 

(4)  conventional  TDR  measurement  to  determine  a  reactive  load  impedance. 

The  homomorphic  deconvolution  method  was  proven  to  be  very  effective  in  separating  convolu- 
tions.    The  method  was  not  limited  by  the  overlap  of  consequent  reflections  in  the  time 
domain. 

In  chapter  4  the  theory  of  the  time  domain  dielectric  spectroscopy  measurements  was 
presented.     Conventionally  this  method  of  measurement  works  only  for  nonoverlapping  reflec- 
tions; i.e.,  reflections  occurring  in  separate  time  domain  windows.     If  the  reflections 
overlap,  only  approximate  solutions  may  be  applied  to  the  frequency  domain  form  of  the 
signal.     By  applying  the  homomorphic  deconvolution  transformation  to  the  TDDS  waveform,  a 
new  method  was  developed  which  was  not  limited  by  overlapping  reflections.     The  new  method 
was  applied  to  the  well  known  Debye  dielectric  model.     This  method  eliminated  the  over- 
lapping reflections  and  determined  the  dielectric  parameters  to  within  0.3%  error  of  their 
model  values.    The  new  method  is  a  general  one  and  has  great  utility  for  time  domain  mea- 
surements on  dielectrics  since  it  does  not  include  any  restrictions  on  the  dielectric 
material  parameters  or  require  the  assumption  of  a  prior  physical  model.     Furthermore,  the 
new  method  is  not  restricted  to  time  domain  spectroscopy  of  only  dielectric  materials  but 
has  utility  for  time  domain  spectroscopy  in  general  such  as  for  magnetic  and  other  phenomena. 
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Chapter  5  began  with  an  introduction  about  the  role  that  the  time  domain  measurement 
techniques  can  play  in  solving  electromagnetic  radiation  problems.     Then  the  distributed 
nature  of  the  expected  cepstrum  in  such  problems  was  discussed  and  compared  with  previously 
studied  discrete  cases.     Two  examples  of  applying  the  homomorphic  deconvolution  transform 
to  an  antenna  system  were  presented,  and  the  transform  proved  to  be  very  successful  in  the 
removal  of  the  unwanted  reflections  and  scattering  components  either  overlapping  or  not 
overlapping  the  main  signal  containing  the  desired  signal  information. 

Chapter  6  was  devoted  to  the  modeling  of  the  feed-through  sampling-head.     The  theory 
and  construction  of  the  sampling  head  were  reviewed.     The  physical  dimensions  of  the  structure 
discontinuities  were  used  to  formulate  an  equivalent  network  for  the  sampling-head  input 
and  output  lines.     The  dimensions  of  the  biconical  cavity  together  with  TDR  measurements 
were  used  to  model  the  cavity  part  of  the  sampler.     The  dc  characteristics  of  the  diode 
together  with  TDR  measurements  vs.  sampling  diode  bias  were  used  to  evaluate  various  diode 
model  parameters.    The  strobing  pulse  was  measured,  and  the  data  corrected  for  both  ampli- 
tude and  risetime.    All  of  these  sampling-head  parameters  were  simulated  in  the  form  of  a 
digital  computer  program  which  was  run  to  obtain  a  theoretical  estimate  of  the  sampling- 
head  step  response.    The  program  was  implemented  in  two  steps  because  of  the  large  number 
of  elements  and  information  required  to  simulate  the  sampling  system.     The  step  response 
was  found  to  vary  with  the  diode  reverse  bias;   the  higher  the  reverse  bias  the  faster  was 
the  step  response,  within  limits.     The  step  response  at  the  manufacturer's  suggested  bias 
and  at  the  bias  for  the  fastest  step  response  was  obtained.     The  homomorphic  deconvolution 
transform  developed  in  chapter  2  was  slightly  modified  in  order  to  be  used  to  deconvolve 
the  sampling-head  response  contribution  from  the  output  of  the  sampling  oscilloscope. 
Deconvolution  of  the  observed  waveform  provides  an  estimate  of  the  actual  (true)  waveform 
at  the  input  port  of  the  oscilloscope. 
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APPENDIX  A 

THE  z-TRANSFORM  AND  THE  DISCRETE  FOURIER  TRANSFORM 

The  z-transform  is  a  generalization  of  the  Fourier  transform  in  the  case  of  discrete 
time  signals  (sequences).     The  z-transform  of  the  sequence  x(n),  where  n  is  the  sequence 
variable,  is  defined  as 


z[x(n)]  =  X(z)  =    I      x(n)  z  ^  (A.l) 
n=-  <» 


where  z[     ]  is  the  z-transform  operator  and  z  is  a  complex  variable. 

The  z-transform  X(z)  does  not  converge  for  all  sequences  x(n)  or  for  all  values  of  z. 
The  region  of  convergence,  i.e.,  the  set  of  zvalues  for  which  X(z)  of  a  given  x(n)  con- 
verges, is  very  important.     It  must  be  defined  in  association  with  the  transform  X(z)  in 
order  for  the  inverse  transform  to  be  uniquely  defined.     Generally  speaking,  the  z-transform, 
eq.  (A.l),  converges  in  an  annular  region  of  the  z-plane  given  by 

0  <_  R_  <   I  z  I   <  R^  £  ~  .  (A. 2) 

For  a  stable  time  domain  sequence,  the  convergence  region  of  the  z-transform  includes  the 
unit  circle  |z|  =  1. 

The  inverse  z-transform  is  given  by 

x(n)  =  2^  ^    X(z)  z""-"-  dz  (A. 3) 

where  C  is  a  counterclockwise  closed  contour  in  the  region  of  convergence  of  X(z)  encircling 
the  origin  of  the  z-plane. 

If  the  complex  variable  z  in  eq.  (A.l)  is  expressed  in  the  polar  form 

z  =  re-J  (A. 4) 

the  z-transform  will  take  the  form 

00 

x(re-^  )  =    I      x(n)  re-'  (A. 5) 

n=-oo 

which  can  be  viewed  as  the  Fourier  transform  of  x(n)  multiplied  by  an  exponential  sequence 
r  If  r  =  1  (|z|  =  1)  the  z-transform  becomes  identical  to  the  Fourier  transform  of  the 

sequence, 

X(ej'^)  =    I      x(n)  e-j^".  (A.6) 
n=-«° 

For  finite  duration  sequences,  with  total  number  of  samples  N,  the  z-transform,  eq. 
(A.l),  will  take  the  form 

N-1  _  - 

X(z)  =    I    x(n)  z  (A.7) 
n=0  "1^^- 


or 

N-1 

X(e"''^)  =    1    x(n)  e  -'"^  ,  on  the  unit  circle.  (A. 8) 

n=0 

For  computational  reasons,  only  a  finite  number  of  samples  of  the  z-transforms  can  be 
calculated.     Sampling  the  z-transform  at  N-equally  spaced  samples  on  the  unit  circle 

|z|  =1  results  in  the  discrete  Fourier  transform  DFT, 

X(k)  =  XCe^"")  ,     =  2-k/N 

,  ,  -i(2TTk/N)n 

=     V    x(n)e  (A. 9) 

n=0 

wnere  k  is  a  sequence  variable  that  takes  on  values  from  0  to  N-1.     By  this  sampling  process 
the  finite  duration  sequence  x(n)  is  forced  into  being  periodic  with  period  N.     Tne  trans- 
form X(k)  will  also  be  periodic  with  the  same  period.     The  discrete  inverse  Fourier  trans- 
form DIFT  is  given  by 

x(n)  =^  Y  X(k)  e^C2.k/N)n_ 
k=0 

Both  tne  z-transform  and  the  DFT  have  the  following  listed  properties  even  though 
these  properties  will  be  presented  only  in  the  z-transform  form.     It  is  to  be  remembered 
tnat  tne  DFT  corresponding  relations  can  be  obtained  by  replacing  the  complex  variable  z  by 
^j(2T:k/N)n  operator  z[     ]  by  DFT[     ].     The  properties  are: 

1.  Linearity 

z[a  x(n)  +  b  y(nj]  =  a  X(z)  +  b  Y(z) .  (A. 11) 

2.  Sequence  shift 

-n 

z[x(n-n^)]  =  z     °  X(z).  ^  (A. 12) 

3.  Convolution 

z[x^(n)   *  x^Cn)]  =  z[x^(n)]    •  z[x2(n)].  (A. 13) 

Finally,  one  additional  but  very  useful  property  of  the  z-transform  will  now  be  discussed; 
it  is  the  multiplication  by  an  exponential  sequence.     This  property  is 

z[a^x(n)]  =  X(z/a).  (A. 14) 

The  importance  of  this  property  lies  in  its  ability  to  scale  the  convergence  region  by 
tne  factor  1/a,  tnus  enabling  the  inclusion  of  the  unit  circle  in  the  convergence  region 
and,  consequently,  stabilizing  unstable  sequences.     This  method  can  also  be  used  to 
remove  the  transform  zeros  that  may  lie  on  the  unit  circle  giving  rise  to  singularities 
upon  taking  the  complex  logarithm.     This  method  of  scaling  the  z-transform  is  often 
called  the  exponential  weighting. 
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APPENDIX  B 
COMPUTER  SUBROUTINES  HDT  AND  HDF 


SUBROUTINE  HDT  (X.H.N. NSH) 
Ctttt  T^C  HOMOMORPHIC  DECOHUOLUTION  TRANSFORM 
C****  FOR  TIME  DOM  SIGNALS 

CmPLEX  X<1).H<1).CLX,CC 
C****  TAKE  T^E  FOURIER  TRANSFORM 

CALL  FFT  (X.N. .TRUE.  > 
C****  CALL  SUBROUTINE  HDF 

CALL  HOF  (X.H.N.NSH) 
C**.t*  INl€RSE  FOURIER  TRANSFORMS  OF  SEPARATED  COMPONENTS 

CALL  FFT  (X.N. .FALSE. ) 

C^  FFT  (H.N.  .F^^SE.  ) 

RETURN 

END 


SUBROUTINE  HDF  (X/H.N.NSH) 
Ctttt  Jm  HCM3M0RPHIC  DECONUOLUTION  TRANSFORM 
C****  FOR  FREQ  DOM  SIGNALS 

COMPLEX  X< 1 >. H( 1 >. CLX. CC 

PI=4.^TAN<1.0) 

N2=N/2-t-l  .  . 

C****  CONTINUOUS  COMPLEX  LOGARITHM 
C=0.0 
ARX=0.0 
00  2©  L=1.N2 

ARS=ARX  - 
CC=X(L) 

cm,  CLG  (CC.CLX) 

ARX=AIM^^CLX> 

AR«*m-^^ 

IF  <m.GJ.'^iy  c=c-2.«Pi 

IF  (AR.LT.-PI)  C=C+2.«PI 

CC3OIPLX<0.0.C) 

X<L>=CLX+CC 
20  CONTIM^ 
Ctttt  LINEAR  PHASE  ELIMINATION 

NSH=AIM^^X(N2))/PI 

TYPE  "    SEQ  SHIFTED". NSH."  SAMPLES" 

00  45  L=1.N2 
45  X(L)»X<L>-CMPLX<0.0.2.*<L-n«NSH*PI/H) 
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K3=N2-2 

DC  25  L=1>K3 
25  X<N2+LM:0NX<X<N2-L>) 
C*m  INUERSE  FOURIER  TRANSFCfRM 

CALL  FFT  <X/N. .FALSE.  > 

TYPE  "  CEPSTRUfI" 
Ctttt  SEPARATION  AhiO  IKUERSE  PROCESSING 

CALL  SEP  CX>H/N.NSH) 

RETURN 

END 


SUBROUTIKE  FFT  (X.N.FI) 
C****  FAST  FOURIER  TRANSFORM  ^  FI=.TRUE. 
C****  INUERSE  TRANSFORM  :  FI=.FALSE. 

CQ^SPLEX  XCN?.U>W>T 

LOGICAL  FI 

EH=N 

M=#:^CG<  EN  >/?LOG<  2 .  >+0 . 5 

N2=N/2 

N1=N-1 

J=l 

DO  7  1=1. Nl 

IF  (I. (2. J)  GO  TO  5 

T=X<  J ) 

>«  J  >=X<  I  ) 

X<  I  )=T 

5  K=N2 

6  IF  <K.GE.J)  GO  TO  7 
J=J-K 

K=K/2 
GO  TO  6 

7  J=J+K 
PI=4.tATAN<l .0) 
DO  20  L=1.M 
L2=2**L 
Ll=L2/2 

IF  <Ll.EQ.i)  GO  TO  8 

U=<i..0.) 

D=L1 

IF  <  .NOT.Fn  D»-0 

W=CriPLX(  COS<PI/D).-SIN<PI/0>> 

8  DO  20  Jal.Ll 
DO  le  I«J.N;L2 
IP1«I-H-1 
T»X<IP1 ) 

IF  <J.^€.1>  T»T*U 

X<  IPl  >=X<  I  )-T 
10  X<  I  )=X<  I  HT 

IF  (Ll.NE.n  U«U*W 
20  CONTINUE 

IF  <FI)  GO  TO  40 

DO  38  I«1.N 
30  XCI>=X<I)/EN 
40  RETURN 

END 
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SUBROUTINE  SEP  <X.H.N.NSH) 
Ctttt  SEP^JRATION  OF  CEPSTRUM  COMPONENTS  AND  INUERSE  PROCESSING 
COMPLEX  X<n.H<l).CC 
DI^ENSION  NS<5).NF<5>.RX<1024) 
LOGICAL  HPF.IS 

00  48  L-l.N 
40  RX<L>=REAL<X<L>> 
Ctttt  LINEAR  FILTER  <  SEPARATION  PROCESS  ) 
30  TYPE  "    READY  FOR  SEPARATION  " 

ACCEPT  "NUf^ER  OF  ELIMINATED  REGIONS  «?  'SN2R 

DO  5  1=1.  N2R 

ACCEPT  "ELIMINATE  POINTS  STARTING  AT  •SNS<I) 
5  <^CEPT  "  AND  Ef^ING  AT  %NF<I> 

NS<NZR+1>=N+1 

HPF=.TRUE. 

NS1«NS< 1 )-l 
90  DO  10  L~1.NS1 

X<L>Cff=>L>«R>«L).0.0> 
10  IF  <HPF)  X<L)«<0.0>0.0) 

DO  16  1=1. NZR 

NSIaNS< I ) 

NFIaNF<I) 

DO  12  LsHSI.NFI 

X<L>=8.0 
12  IF  X<L)«RX<L>-X<L) 

NFIls^dHl 

NSI1«NS< I+l >-l 

DO  14  LaNFIl/NSn 

X<L>«aiPLX<RX<L).0.0) 

14  IF  <hPf>  X<L)«»<0. 0.0.0) 

15  CONTINl€ 

OPm  IWERSE  PROCESSING  '  FOURIER  TRANSFORM 

CMl.  FFT  <X.N.  .TRIE.  > 
Ctttt  RETAIN  LINEAR  PHASE  AND  COMPLEX  EXPONENTIAL 

00  20  L«1>N 

CC=CfW«  0 . 0. 2 .  «<  L-1  )*NSH*PI/N  ) 

IF  (  NOT.^PF)  X<L)«X<L)+CC 

X<L>«CEXP<X<L)) 
20  COHTIKE 

IF  <  NOT.HPF)  GO  TO  99 
Om*  HPF  F^Q  DOM  SIGNAL 

DO  70  L^l.N 
70  H<  L>=X<H  > 
Ontt  REPEAT  FOR  LPF  COMPONENT 

HPF=.Fj^SE. 

GO  TO  9® 

99  ACCEPT  "sWTHER  SEP  TRIAL  ?  -laYES.  0-NO  MS 
IF  <IS)  GO  TO  30 
RETURN 
END 
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APPENDIX  C 


COMPUTER  PROGRAMS  TDNA  AND  SMPL 

C.l    Program  TDNA 


C**:*:*  PROGRAM  TDNA  '  TI^E  OOM  NETWORK  ANALYSIS 

C****  THIS  PROGRAM  CAN  HAmLE  ANY  OF  THE  FOLLOWING  NETWORK  ELEMENTS  •■ 
Ctttt    1=  CUR  SRCE..  R=  RES.  0=  DIODE  <NLR).  C^  CAP.  N'  NL  CAP. 
C****    V-  ULTG  SRCE.  L=   M.  T=  TRANSM  LINE 

DI^ENSION  CK651.  n  >.C1T<  11  >.C2<651. 11  ).C2T<  11  >.R0<  11  ).NT<  in 

DIMENSION  IBNCH<  7. 100  >. ELUALC 100 ). NOOOUTC  2. 30  ). UOUT<  30 ) 

DIMENSION  TIME<651>.TDU<651.1).TDUN<1) 

DIMENSION  G<35.35).C<35).U<35>.USAU<35).NU<10> 

DIMENSION  UPLOK  651 . 26  > 

COWON  mMPTS.TMIN.TMAX 

NDIM=35 

M=0 

C^^M^^    READ  If"^UT  DATA 

OUJ-  RB^IBNCH.ELUAL.^CL.NN.NODOUT.NOUT.R0.LL.NU) 
T=TI*ttX<'i'^Jf^TS 

NU1=NU<1)  $  NOUTNU«NOUT-»+IUl  $  NUPl-NUl+l 
LLL=2^  $  NOUTLL«NOUT+LLL 
NUMPTS#^UMPTS-H 
DO  2  Jal.NUMPTS 
TDU<J.1)=0.2S 
2  TIME<J)=<J-1)*T 
TDU<1.1)=0.0 
IS=1 

CALL  PLOT  (TIME.TOtXl.D.NUMPTS.ieWINPUT  US  T) 
Ctttt    ADJUST  ELEMENT  LIST 

CALL  PREPAR<IBNCH.£LUAL.NEL/m.T.NT.R0> 
DO  18  J«1.NN 
USAU<J)a0.0 
18  U<J)»0.0 
DO  2©  L«1.LL 
C1T<L>«0.0 
20  C2T<L>e.0 

WRITE  <6.89)  NEL/W.U-.NWTS.Nyi.NCSJT/NOUTW/TOfrLL.NUPl/LLL 
l«ITE  <6.92>  <<<imCH<JJ.n>.JJ»1.7>.ELUAL<n>)/n-l/Na.> 
DO  50  N«i.NUMPTS 
cm*    SET  UP  NODE  EQUATIONS 

DO  22  n«i.is 

22  TDUH<n>=TOU<N.  ID 

24  CALL  NOOEQ  < IBNCH.ELUAL.NEL.G.C.U.NN.NDIM.USAU.T1ME<N)/T.TDUN/ 
X  C1T/C2T) 
IF  <N.GT.1.0R.M.GT.l)  GO  TO  28 

tRITE  <6.92>  <<<IBNCH<JJ.II).JJ«1.7>,ELUAL<II)>.II-1.NEL) 

DO  26  11=1. NN 
26  WRITE  <6.91>  <G<  II.  JJ).  JsJ=l.NN).C<  II  ).U<  II  ) 
28  CONTINUE 
C«m    SOLUE  NODE  EQUATIONS 

CALL  SLNEQ  <NN.NDIM.G.C) 
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Ctt^t    y<J)-C<J>  IS  THE  INCREMENT  IN  THE  JTH  NODE  UOLTAGE 
ERR=0.® 
DO  30  J=1.NN 
ERR=€RR-K  U<  J  >-C<  J  >  )*3|c2 
30  U<J)=€<J> 

Cmt    TEST  COHUERGENCE  CRITERION 

IF<ERR.GT.l.E-6.AND.ri.LT.15)G0  TO  24 
M?ITE  <6.93)  N.U<8> 
CALL  OUTPUT<NODOUT,NOUTLL/U.UOUT> 
00  32  IIsl/NOUT 
32  UPLOT<N.II)«U0UT<n) 
DO  34  JJ«!/NU1 

n»wuT*jj  $  NUjaNU<jj+n 

34  UPLOTCN/IIMKNUJ) 

TO  36  Lal/LLL 

II»NOUTNV+L  $  JJ-NOUT-H. 
36  UPLOT<N.II>»UOUT<JJ) 
Cmt    STORE  TRANSMISSION  LINE  CURRENTS 

00  38  L«1/LL 

IL«N0UTNU+2*L-I 

C1<N.L><^PL0T<N.  IL>^^L>-C2T<L> 
38  C2<N.L>«UPL0T<N.  IL+1  VR0<L)-C1T<L) 
Ctttt    STORE  NODE  UOLTAGES  FOR  UPDATING  CAPACITORS 

DO  40  Jsl.NN 
40  USAU<J)=U<J> 
M=0 

DO  45  L=1.LL 

IF  <N.LT.NT<L>)  GO  TO  50 
Ontt    UPDATE  TL  DELAYED  CURRENT  SOURCES 

^«^i=N~HT<LHl  S  IL«N0UTW'»-2$L'-1 

C1T<  L )=UPLOT<  NNTl . IL )/R0<  L )+Cl<  NNTl . L ) 
45  C2T<L>=W-0T< NNTl.  IL+1)/1^LHC2< NNTl /L) 
50  C43NTINIE 

00  55  II=1.N0UTNU 

ITITL=IT<n.NOUT> 
55  m±  PLOT  <TIfCUPL0T<l.II),Mf1PTS.ITITL) 

I^ITE  <6.94>  TMAX 

IF  <NOUT.EQ.0>  GO  TO  65 

WRITE  <6.95)  <JJ.<N0D0UT<II/JJ)/n«1.2>/JJ«l.N0UT> 

DO  60  Nal.NUMPTS 

TNR=TIPE<NVTMAX 
68  WITH  <6.91)  TNR.<UPLOT<N/II>.n«l.NOUT> 
65  IF  <NV!.EQ.0)  GO  TO  80 

IF  <NOUT.hE.0>  WRITE  <6.94>  Jfm 

WRITE  (6/96)  <NU<  ID.  n«2.NUPn 

DO  70  N-l.NUMPTS 

TNRaTIMECNVTMAX 
70  WRITE  <6.91)  Th«.<UPLOT<N.n-^HaiT).II"l/NUl) 
80  WRITE  <7.97)  <UPL0T<N.N0UT+1 >.N»iiNU«PTS) 

89  FORMAT  (1016) 

90  FORMAT  (8F10.0> 

91  FORMAT  <1X/10F13.7) 

92  FORMAT  < IX. A4/6I3/E15.5> 

93  FORMAT  <I7.E1S.9) 

94  FORMAT  <1H1.«  Tm<  »  «/E12.5/> 

95  FORMAT  <8X.«T/0rMAX    «.8<«     BW/ 12.»  0ET«.2I2>> 

96  F^WAT  <^  T/TMAX    *.9<*  N0DE»/I3)) 

97  FORMAT  <5E16.9> 
STOP 

END 


c.2    Program  SMPL 


Ctttt  PROGRAM  SMPL  =  SAMPLING-fCAO  STEP  RESPONSE 

DIMENSION  IBNCH<7.42).ELUAL<42) 
DIMENSION  G<  16. 17  ).  C<  16    U<  16  >. USW<  IS  > 
DIMENSION  TDU<500.2>.TDUN<3).TOC»S<500> 
Dlf^SION  CK 500 >.C2< 500 )..UPLOT< 500.2) 
M=0 

Cttt^    READ  INPUT  DATA 

TYPE  "INPUT  ■  NEL.NN..NPTS,NUf^TS.NT/T.R0.NBB.Nll.N2i>N12.N22" 
ACCEPT  N£L.NN.NPTS.NUMPTS.NT>T.R0.N^,Nll.N21.Ni2.N22 
ACCEPT  "INPUT  =  NBC^..NC..U0.T0UN<3)  ".NBGN.NC.U0.TDUN<3) 
READ  (13.92)  <<  IBNCH<  II.  JJ>..  II=1.7).ELUAL':  JJ).  JJ-l.NEL) 
WRITE  <10.92)  <<IBNCH<II.JJ).II=1.7).ELU^<JJ).JJ=1.NEL) 
READ  <  13.97)  <TDU< J..2). J=1.NPTS) 
READ  <  13.97)  <TIXK J.  1  ). J=1.NUMPTS) 
TIXi2F=TDU<NPTS.2) 
DO  4  J=i.NUMPTS 
TDU<  J. 1 )=TDU<  J. 1 )/25 . 0 
IF  <J.GT.I^TS)  TDU<J.2>«TDU2F 
TDU<  J.2)=TDU<  J.2)/R@ 
4  TIXJS<  J)=TDU<  J.2) 
D1IU=£LUAL<NBB) 
D2IU=ELUAL<NB8+1) 

ACCEPT  "FEED  SACK  FACTOR  =  ".F^ 

NSTP=i 

NPTS=250 

DO  80  NSPE=NBGH.Wjr^TS.NSTP 

W3«NSPE-I^TS 

IF  <ND.£Q.0)  GO  TO  16 
8  J«1.ND 
8  TDV<J.2)«TDUS<1) 

DO  10  Jal.NPTS 
10  TDU<J+ND>2>«TDUS<J> 
16  CONTINUE 

DO  18  J«1.NN 

USAU< J>=0.0 
18  U<J)=6.9 

ELU!^<       )Oi  lU-^TD*^  3  )*1 .  E5 

ELUAL<  hm-i-i  >=02IU-TDW<  3  )*1 .  ES 

ciT»0.a 

C2T«0.0 
DO  50  N=i.NSPE 
C****    SET  UP  NODE  EQUATIONS 
TDUN<l>»TDU<N.l) 
TIX^2)=TDU<N.2) 
TIME=<N-i)*T 

24  CALL  NOKQ  <  IBNCH.ELUAL.NEL.G.C.U.NN.NDIM.USA<J.T2f€<N)>T.TDUN. 
X  C1T.C2T) 
CW**    SOLUE  NCOE  EQUATIONS 

CALL  SLHEQ  <NN.NOIM.G.C> 
cm*    U<J>-C<J)  IS  THE  INCREMENT  IN  THE  JTH  NODE  UOLT<^ 
ERO=0.0 
DO  30  ^l.m 
ERO=€RO+<  U<  J  >-C<  J  )  >**2 
30  IXJ>=C<J> 
f1=fHl 
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amt    TEST  CONUERGENCE  CRITERION 

IF<ER0.GT.1.E-6.AND.M.LT.15)G0  TO  24 

W.OT< N.  1  )=U< Nl i  )-U< N21  > 

UPLOK  N/  2  >=U<  N12>-U<  N22  ) 
CUtt    STORE  TRANSMISSIC»^  LI^€  CURRENTS 

Ci<  N  )sUPLOT<  N. 1 )/R8-C2T 

C2( N  >=UPLOT<  H. 2 )/R0-ClT 
CUtt    STORE  NODE  UOLTAGES  FOR  UPO#^TING  CAPACITORS 

DO  48  J^i/NN 
40  USAU<J)sU<j) 

IF  <N.GT.l)  GO  TO  42 
UCiIaU<2)HKl) 
UC2I»U<a>-U<9) 
42  Ma® 

IF  (N.LT.NT)  GO  TO  50 
Ctttt    UPDATE  TL  DELAYED  CURRENT  SOURCES 
NNTl^N'-NT'^l 

C1T=UPL0T<NNT1. 1 )/R0+Cl<NNTl > 
C2T=UPL0T<  NNT 1 . 2  )/TOK:2<  NNT  1  > 
50  CONTII'^^ 

USC=<  U<  8  MX  9  )  )-UC2I-<  U<  1  )-U<  2  )  )-UCl  I 
NC=hC+l 

TOUN<  3  )=TDUN<  3  )+FBF:*sUSC 
NCT=<  NC-1 >*T/1 . E-12+ . 5 
U0=U0+USC 

mUE  <10..93)  NC.NCT.U0/TaiN<3) 
80  CONTINUE 

92  FC^MAT  <1X.A4.I6.5I3.E15.5) 

93  FORflAT  <1X.I5.5X.I5.2X>"PS".2E17.6) 
97  FORMAT  <5E16.9) 

STOP 
Eh© 
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C.3     Subroutines  Used  in  TDXA  and  S>fPL 


SLBROUTINE  RE^^<  I^O^.  ELUAL.  ^CL.  NN.  NODOUT.  HOUT.  R^.  LL.  ) 
DIMENSION  IBNCH;?.  1  ^^ELUALC  1  >.NOD0iIT<2>  1  )>DUmV<  15).Re<  1  > 
INTEGER  TA.NUC  1  ).  IU25).NUT(  10);NUA<9) 
COMMON  NUt1PTS;TFMIH.TFmX 

DATA  IBLAN*</1H  ^  ^ NUi^^/lHA.  IHe.  IHC.  IKO.  1H£.  IHF.  IHG/ 1HH>  IHI/ 
C****    REM)  AND  WRITE  NAUE  CARD 

READ  <5..50)  (OmMYC  J).  J=i.  15? 

l^ITE  (6.52)  <DUmY(  J).  J=l.  15) 

L=8  $  NN=e  *  N£L=5  $  NCHJT=0 
Ctttt    R£f€  hCXT  ELEMENT 

READ  (5.54)  <  IBNCH<  J. NEL).  J=l .  5 ).  ISGN;  lOUT.  IS.  < DLWYC  J ).  J=l . 2). 
X  <N<vT(  J).  J=l.  10) 

C****    READ  FIRST  CHARACTER  CF  NAf^  <=ELEr£NT  TYPE) 

BACKSPACE  5 

READ  <5.5€)  TA 
OPm    DECODE  ELEflENT  TYPE 

NTYPE=NOX€<  TA  ) 

IF  <NTYP£-9)  20.35. 10 
Ctttt    SAC€  GEN^ERAL  INFORMATION  CARD  DATA 
10  N!J1PTS=IBf^H<4.N£L) 

TFMAX=DLfim'<  1  )  $  TFNIN«CX«MY<2) 

IBNCH<  l.h€L)=4H 

NEL=NEL-1 

DO  15  J=1.10 
15  NU<J)=NJT<J) 

NU1=NU<1) 

IF  <NUl.e3.0)  GO  TO  5 
fcRITE  <6.58)  <NivKJ+l).J-l.NUl) 
WRITE  (€.60)  J).J«1.NU1) 
GO  TO  5 
20  ELUALCNELXBJ««Y<1) 
Opm    STORE  TYPE  OF  ELE^E^f^ 
IBNCH<6.Ha-)«NTYPE 
IBNCH<7.S£L)==0,0 

IF  (CNTYPE.EQ. 1 ,0R.NTYPE.EQ.6).AN0. IS.NE.0>  IBNCH< hSL >«IS*4 
cm*    STORE  OUTPUT  NODE  NUMSERS 

I0UT=NOXE<  I  OUT) 

IF  <I0UT.^€.6)  GO  TO  30 

NOUT=NOUT+1 
cm*    CHECK  SIGH  OF  CBJTPUT  UOLTAGS 

IF  <I9>I.E]3.IBLANK)  GO  TO  25 

NOOOUT<  1 . NOUT >=  iaCH<  5 >  hCL) 

NOCXXIT<  2.  NOUT  >=IBNCH<  4.  NED 

GO  TO  30 
25  HOOOUTC  l.N0UT>*I^CH<4.hCL) 

NOOOUTC 2. NOUT )« I BNCH< 5 . NEL ) 
cm*    STORE  HIGHEST  NODE  NUMBER 

30  NN=*W<9<NN.IBNCH<4.fCL).IBNCH(5.h€L)) 
cm*    READ  OUTPUT  NODES  OF  TRANStllSSION  LINE 

IF  <IBNCH<6.NEL).NE.8)  GO  TO  5 

L=L+1 

READ  <5.62>  N3.N4.R8<L) 
hfW1AX0<hM.N3.N4) 
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IBNCH<2.NEL)=N3  $  IBNCH<3.he.>«N4 

GO  TO  5 
35  IBNCH<l/hEL)«4H 
N£L=NEL-i 
LL=L 

DO  40  Lsl^LL 

liL=ll<L)  $  IL»NOUT+2:5cL-i 

N00ajKl/IL>«IBNCH<4>nL>  $  NOOOUT<  1.  IL-H  )«IBNCH<2.  IIL) 
N000UT<2.2L)«IBNCH<5.nL>  $  N0D0UT<2> 2L+1  >«IBNCH<3. UL> 
40  CONTINUE 
URITE  <€/64) 
L«0 

QQ  4^  Jsa  I  ,  NEL 

URITE  <6>€6>  ieNCH<i/J>.IBNCH(4/J>/IBNCH<5.J>^ELMAL<J>/ 
X  2BNCH<3.J)/IBNCH<2/J> 

IF  <IBNCH<6.J).EQ.8>  WRITE  <6.68)  IBNCH<2.  J).  IBNCH<3.  J).R0<L-#-l  ) 
45  CONTINUE 

RETURN 
50  FC^^IAT  <  15^) 
52  FORMAT  <  10<  IH*).  iS^H.  !0<  1H«)) 

54  FORMAT  <A4. A1.A4. 14. 13. 1X/A1.A1> I1.2E15.5/ 10I3> 
56  FORMAT  <A1 > 

58  FORMAT  OUTPUT  NODES  AND  THEIR  CORRESPONDING  ALPHABETIC 

X  *NA^€S  USED  IN  MICROFILM  PLOTS*/. 1X.9<2X. I3>> 

60  FORMAT  <iX.9<4X.Al)) 
62  FCm^T  <10X.2I3.3X.E15.5> 

64  FCFMAT  <//*  BRANCH  h03E       ELEMENT       BRANCH  CON-*/ 
X  *  NAME     NOS.        UALUE  NAME  TROL«) 

66  FORMAT  < IX. A4.2I3.E15.5. IX. A4.3X. Al > 
68  FORMAT  <5X.2I3.E15.5) 
END 


SUBROUTIfC  PREPAR<IBNCH.ELUAL.NEL.NN.T.NT.R0> 
DIMENSION  IBNCH<7.1>.ELUAL<1).NT<1).R0<1> 
COMMON  K01PTS.TMIN.TMAX 
DATA  IBLANK/IH  /. IREF/lHI/ 
C*m    mTCH  CCOT50LLED  Pm  CONTROLLING  ELEMENTS 
00  2  I«|.NEL 

IF  <IBNCH<6.I).NE.1>  GO  TO  2 
IF  <IBNCH<2.I>.EQ.IBLANK>  GO  TO  2 
IF  <IBNCH<2.I).EQ.IREF)  GO  TO  2 
DO  1  J^l.NEL 

IF  <IBNCH<3.I>.NE.IBNCH<1.J))  GO  TO  1 
cm*    STORE  CONTROLLING  NODE  NUMBERS  FOR  USE  IN  NOOEQ 
IBNCH< 2. I )=IBNCH< 4. J ) 
IBNCH< 3. I >«IBNCH< 5. J > 
GO  TO  2 

1  CONTINUE 

2  CONTINUE 
f€Ll=NEL 

C**«*    MODEL  INDEPENDENT  UOLTAGE  SOURCE  BY 
C  GVRATOR  AND  INDEPE^«DENT  CURRENT  SOURCE 

00  3  I»l/NELi 

IF  <I8NCH<6.I).NE.S>  GO  TO  3 
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IF  <iaCH<2.I).NE.IBUW)  GO  TO  3 

GYRATR<I8HCH.ELU^.N£L>NN.M.0.1.8) 
cm*    MO»v€        CHfii^  TYPE  OF  IbCEFEHDEHT  SOJRCE 
IBNCH<4.I)=NN 

3  CONTINiJ€ 

Om*    f-QOEL  INDUCTOR  BY  GYRATOR  ^  CAPACITOR 
DO  5  I-l.NELl 

IF  <ISHCH<S.  n,NE.7>  GO  TO  5 
C^^  GYRATR<ISNCH.ELUAL.hCL.NN/M.@;l.e) 
cm*    flOUS  A^C  REPLACE  INDUCTOR  WITH  CAPACITOR 
IBNCH<4.I>«NN 
IBNCH<5.I)»e 
IBNCH<6. I )=4 

5  CONTINUE 

C*»*    MODEL  SHORT  CIRCUIT  BY  A  GYRATOR  Afd  OPEN  CIRCUIT 
DO  6  >l.hCLl 

IF  <IS^>K6.I).NE.2)  GO  TO  6 
IF  <EL<^CI).NE.0.0)  GO  TO  6 

SCCT  <IBNCH.ELUAL.h€L.h*^;I) 

6  COfTINUE 

C*«*    REPLACE  TRANSMISSION  LINE  BY  ITS  MODEL 
L=8 

1X3  7  I=l.hELl 

IF  <I&CH<6.I).h€.8)  GO  TO  ? 
L=L-»-l 

NT<L)=£LUAL<I  VT+e,5 

CALL  TRLIN  < IBNCH. ELUAL. NEL> I . Re< L ) ) 

7  CONTINLE 

Ctttt    MODEL  NONLINEAR  CAPACITOR  BY  ITS  COMPANION  MODEL 
DO  8  I«i.NELl 

IF  <IBNCH<6.I),hE.5)  GO  TO  8 
Tl=l./T 

CALL  GYRATRC IBNCH. ELUAL. NEL. NN/ I . Tl . 1 . 0 ) 

IBNCH<7.hEL-l >=1 
cm*    CHAM3E  TO  UARIABLE  CIRRENT  SOURCE 
C  NOTICE  SIGN  CHANGE  IN  hODEQ 

IBNCH<3.I)«NW 

IBNCH<6.I>«1 

I»CH<7.I)»1 

C****    ADD  UARIABLE  COhOUCTANCE  AND  CURRENT  SOURCES 
NEL«NEL+1 

IBNCH< 1 . NEL >«IBNCH< 1 . I > 

IBNCH<2/h€L>«lH 

IBNCH<3.NEL>-I 

IBNCH<4.h£L>«NN 

IBNCH<5/NEL>»e 

I8NCH<6.hCL>«2 

IBNCH<7.^€L)=1 

NEL-h€L+l 

IBNCH< 1 . NEL >«IBNCH< 1 / n 

IBNCH<2.h€L>»lH 

IBNCH<3.NEL>»NN 

IBNCH<4.NEL>=NN 

IBNCH<5.NEL>=0 

IBNCH<6.NEL>=1 

IBNCH<7.NEL>=2 

8  CONTINUE 
RETURN 


SUBROUTIfC  NOOEQ<IBNCH.ELUAL.^CL.G.C.U.^#^>NOIM.USW.TI^fE.T.TDUAL. 
X  CiT.C2T) 
REAL  LAMDA 

DI^E^^SION  IBNCHK  7.  l  >,ELUAL<  1  ).TOUAL<  1  ).C1T<  1  ).C2T<  1 ) 
DIMENSION  G<N0II1. 1  ),C<  1  ).U<  1  ).USAU<  1) 
COmON  NUMPTS..TI1IH.Tf1AX 
DATA  IK-ANK/IH      INIT/9// IREF/lHI/ 
DATA  CSAT/'2.76E-ll/'.LAMDA/'38.6i/' 
C^^c***    ZERO  WOE  mTRIX  Am  CURRENT  UECTOR 
DO  5  K=l.m 
C<K)=0.0 

DO  5  j=i.m'  • 

5  G<J..K>=5@.8 
L=0 

00  80  J==1.NEL 
Om:t    Nl  AM)  N2  ARE  NODE  WJ-SERS 

N1=IBNCH<4.J> 

N2=IBNCH<5.J) 

BU=e.0  $  BUN=0.0 

IF  <N1.EQ.0)  GO  TO  10 

B«v>=U<Ni>  $  BUN=USAU<N1) 
10  IF  <N2.£Q.0)  GO  TO  15 

BU=BU-U<N2>  $  BUNaBW-USAU<N2) 
15  CONTINUE 
CtUt    DETERMINE  ELEHENT  TYPE 

NTYPE=IBNCH<6.J> 

GO  TO  <2®.25/70.7S>.NTYPE 
Ctm    TEST  TO  SEE  IF  INDEPENDENT  CURRENT  SOURCE 

20  IF  <I»CH<2.J).EQ.IBUW)  GO  TO  35 
C*t««    MODIFY  NONLINEAR  CAPACITOR  MODEL  AT  DC 

IF  <TIME.EQ.0.0.AND.IBNCH<7.J).EQ.l)  GO  TO  8® 
Omt    LOAD  UOLTAGE  CONTROLLED  CURRENT  SOURCE  INTO  G 

CfiLL  UCCS<Nl.N2.IBNCH<2^J)/IBNCH(3>J>iELUAL<J>^G>N0IM> 

GO  TO  S© 
Omt    LtM>  CQf«XTANCE  INTO  G 

25  IF  <IBNCH<7.J).EQ.0>  GO  TO  30 

ICNL«IBNCH<3>J> 

NC1=«IBNCH<4.ICNL>  $  NC2«IBNCH<3.  ICNL> 
8MC=0.0 

IF  <NC1.ME.0>  BUC«U<NC1> 

IF  <NC2.f^.0>  BUC=BUC-U<NC2> 

ELUAL<J>=CNL<BUC) 
30  COM5=i.0/€LUAL<J) 

CALL  BMCH<N1.N2.C0ND.G>NDIM) 

GO  TO  80 
35  IF  <IBHCH<7.J).GT.4>  GO  TO  60 

ICNT=ie^CH<3.J> 

GO  TO  <  65. 40. 45  . 50. 55    IBNCfK  7/ J  HI 
cm*    MCS3IFY  NLC  MODEL  AT  DC 

40  IF  <TIME.EQ.0.0)  GO  TO  8® 
C*m    CmHGE.  SIGN  AS  MENTIONED  IN  PREPAR 
ELUAL<J>«-VSAU<ICNT)/T  $  GO  TO  63 
45  EL«M.<J>aeUC-CONO«U<ICNT)  $  GO  TO  65 
58  L=L*1  $  ELUAL<J)«C2T<L)  $  GO  TO  65 
55  ELUAL<J>sCiT<L>  $  GO  TO  65 
60  IS«IBNCH<7.J>-4  $  ELUAL<  J)»TDUAL<  IS) 
65  CALL  SRC£<N1.N2/ELUAL<J).C) 
GO  TO  ^ 
Ct$$t    DEFINE  DICDE  PARAMETERS 
cm*    ^m.Y  INITIAL  DIODE  IF  ON  FIRST  ITERATION 
70  IF  <INIT.EQ.0>  BU»€LUAL<J) 

-194- 


C3*»    IMPOSE  HEURISTIC  LIMITATION  OF  .8  UOLTS  ON  DIODE  BU 

BU=A«IN1<BU.0.8) 

DEX»0.0  S  BUJl}«BU*LArClA 

IF  <SULMD.GT.-6?S,@)  C€X»EXP< BULMD ) 

GO«CS^T«IJVtM«DEX 

CUR-CSAT*<  OEX-1 .  @  >-iU*0 
cm*    LOAD  (2)  AND  CUR  INTO  G  C 

C<^  BHC>4<N1.N2/GD/G.N0IM) 

CALL  ^CE<N1.N2.CUR.C) 

GO  TO  88 

Ctttt    FORM  CJ^AHIOH  flOOEL  OF  CAPACITOR;  SKIP  IF  DC 
75  IF  <TI?«.EQ.8,8)  GO  TO  88 
GC«ELU^^<  J  )/T 
CUR=8UN^ 

CALL  ^40KNl.N2.GC.G.NOIf1) 
CALL  SRC£<N2.N1.CLR.C) 
80  CONTIhJLS 
INIT=1 
f^TURN 
END 


aJ8R(XJTI^€  S_NEQ<N>M3IM.A.B) 
Ctttt    SOLUTION  OF  H  LIh€AR  EQUATICNS  A< I . J )*B< J >«B< I > 
DIMENSION  A<NDIM.1)>B<1) 
N1=N+1 
DO  18  1=1. N 
10  A<I.N1>=B<I> 
DO  40  K=1.N 
K1=K+1 
G=8.8 

DO  20  I=K.N 
P=ABS<A<I>K)) 
IF  (P.LE.G)  GO  TO  20 
G=P  $  L=I 
20  CONTINL€ 
DO  25  J=K. Nl 

G=A<L>J)  $  A<L. J>«A<f<. J)  $  A<K.J« 
25  CONTIWJE 

P=A<K.K) 

DO  38  J=K1.N1 
38  A<K.  J)=A<K.  J  VP 

DO  48  I»1.N 

IF  CI.EQ.K)  GO  TO  48 

P=-A<I.K) 

DO  35  JaKl.Nl 
35  A<I.J)«A<I.JHP*A<K.J) 
48  CONTINUE 

DO  45  >1/N 
45  B<J>«A<J.N1) 

I^TLRN 

END 
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SUBROUTIhE  GYRATR< IBNCH.ELUAL.NEL.NN, I.G1/G2> 

DIMENSION  ELUAL< 1 >/ IBNCH<  7.1) 
Ctttt    THIS  SUBROUTINE  GENERATES  TWO  UCCS'S  TO  SI^aJLATE  GYRATOR 

N1=IBNCH<4.I)  *  N2*IBNCH<5/n 

NN=NN+i 
Ctttt    ADD  FIRST  UCCS 

IBCH<  1 .  ^€L  >«IB^CH<  1 ..  I  > 
IBNCH<2.hEL>=NN 
IBNCH<3.hEL>=0 
IBNCH<4.NEL)=N1 
IB^CH<5.hEL)=N2 
IBNCH<6.NEL)=1 
IBNCH<?.NEL>=€ 
ELUAL<NEL>=G1 
Ctttt    ADD  SECOND  UCCS 
^CL=NEL+l 

iaCH(  1  > f€L  )=IBNCH<  M  > 

IBNCH<2.NEL)=N1  . 

IBNCH<3/NEL)=N2 

IBNCH<4.NEL)=0 

IBNCH<5.NEL)=NN 

IBNCH<€.NEL>=1 

IBHCH<?.NEL)«0 

ELUAL<NEL)=G2 

RETURN 

END 


SUBROUTINE  SCCT< IBNCH.ELUAL.NEL.NN, I) 
DIMENSION  ELUAL< 1 )/ IBNCH<  7.1) 
Ctttt  THIS  SUBROUTINE  REPLACES  SC  BY  GYRATOR  AND  X 
Ni=IBNCH<4.I)  $  N2«IBNCH<5/I) 
NN=NN+1 
IBNCH<2.  I>=NN 
IBNCH<3.I>=0 
IBNCH<6.I>sl 
IBfrCH<7.I>s8 
ELUAL<I>=1.® 
fCLsNEL+l 

I8NCH< 1 , NEL )«IBNCH< 1 . I ) 

I^CH<2.hEL>=Nl 

IBMCH<3.hEL>=N2 

I0NCM<5>ICL)«NN 

IKCH<6.f€L)=l 

I^CH<7,NEL)a0 

ELUAL<NEL)«S.0 

RETURN 

END 
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SUBROUTINE  TRLINC IBNCH.ELUAL/NEL. I.R0> 
DIMENSION  ELUAL< 1 ). IBNCH<7; 1 ) 
N1=IBNCH<4. 1  )  *  N2-IBNC-HC5;  I  > 
N3=IBHCH<2. 1  )  $  N4=I^H<3.I) 

is^b:h<2.  I  >=1H 

I»CH<3.I)=iH 
IBNOK6. 1  >=2 
IBNCH<7.I>=0 
ELiv¥^<  I  >=R0 
^CL=^€L+l 

IBNCH<  1 . ^£L >=IBNCH<  1 . 1) 

IBNCH<2.h€L>=lH 

IBNCH<3.^€L>=1H 

IBfCH<4.h€L>=N2 

iaCH<5.hEL)«Nl 

IBNCH<6.^CL>=l 

IBNCH<7.NEL>»3 

ELUAL<KEL)=1.8 

NEL=NEL+1 

IBNCH<  1 . hEL)=IBNCH<  1 . 1 ) 

IBNCH<2.NEL)=1H 

IBNCH<3.NEL>«1H 

IBNCH<4.NEL>«N4 

IBNCH<5.NEL>«N3 

IBNCH<6/NEL>«1 

IBNCH<7.NEL>«4 

ELUAL<NEL>ii  .8 

fe.=NEL+l 

ISNCHC  1  > NEL)«IBNCH<  1 . 1 ) 

IBNCH(2/NEL>«1H 

IBNCH<3.f€L>»lH 

IBNCH<4.NEL>«N3 

IBNCH<5.NEL>^ 

IBNCH<6.h£L>2 

IBNCH<7.NEL>=0 

ELUAL<NEL)=R0 

RETURN 

END 


SUBROUTINE  UCCS< N1.N2>N3.N4.C0N0>G;NDIM) 
DIfCNSIOH  G<NOIM.n 

IF  <Nl.h£.0.AND.N3.NE.3)  G<Nl/N3)«G<Nl.N3)-*0>© 

IF  CN2.Ne.0.AND.N3,h€,0)  G<N2.N3>«G<N2.N3)-C0fCl 

IF  <Nl.NE.0.^.N4.hE.0)  G<Nl.N4)«G<Nl.N4)-C0^Cl 

IF  <N2.NE.0.f^.N4.NE.0)  G<N2.N4>«G<N2.N4>'H3>C 
RETURN 
END 
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SUBROUTIhC  BNCH<N1.N2,C0ND.G.NDIM) 

OlflENSION  G<NOIM..  1) 

IF  <Ni.h€.e>  G<N1.N1)=G<N1/N1>-»C0N0 

IF  <N2.HE.0>  G<N2..N2>=G<N2.N2HC0ND 

IF  <Ni.EQ.0.OR.N2.EQ.0)  RETURN 

G<N1.H2>=G<N1,N2)-C0ND  $  G<N2.N1  >=G<N2.N1  )-aMJ 

Em 


SUBROUTINE  SRCE<N1.N2. SOURCE. C) 
DIMENSION  C<1) 

IF  <Nl.r€.0)  C<N1)«C<N1)-S0URCE 
IF  <N2.HE.e)  C<N2>=C<N2HS0URCE 


FU^CTIW  NCODE<NAriE) 
DIMENSION  IREFC9) 

DATA  IREF/lHI . IHR. IHO. INC. IHN. IHU. IHL. IHT. IKE/ 

DO  1  1=1.9 
hCODE^I 

IF  <N^«C.EQ.IREF<n)  GO  TO  2 

1  CONTINUE 
HCODE=10 

2  RETURN 

Em 
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